Merge pull request #2362 from KaiSzuttor/lb_thermostat

LB thermostat test

testsuite/python/CMakeLists.txt

@@ -166,6 +166,7 @@ python_test(FILE field_test.py MAX_NUM_PROC 1)
 python_test(FILE lb_boundary.py MAX_NUM_PROC 2)
 python_test(FILE lb_streaming.py MAX_NUM_PROC 4)
 python_test(FILE lb_shear.py MAX_NUM_PROC 2)
+python_test(FILE lb_thermostat.py MAX_NUM_PROC 2)
 
 if(PY_H5PY)
   python_test(FILE h5md.py MAX_NUM_PROC 2)

testsuite/python/dpd.py

@@ -18,12 +18,13 @@
 #
 # Tests particle property setters/getters
 from __future__ import print_function
-import unittest as ut
-import espressomd
 import numpy as np
-from time import time
+import unittest as ut
 from itertools import product
+from time import time
 
+import espressomd
+from tests_common import single_component_maxwell
 
 @ut.skipIf(not espressomd.has_features("DPD"), "Skipped because feature is disabled")
 class DPDThermostat(ut.TestCase):
@@ -43,11 +44,6 @@ def tearDown(self):
         s = self.s
         s.part.clear()
 
-    def single_component_maxwell(self, x1, x2, kT):
-        """Integrate the probability density from x1 to x2 using the trapez rule"""
-        x = np.linspace(x1, x2, 1000)
-        return np.trapz(np.exp(-x**2 / (2.*kT)), x)/np.sqrt(2.*np.pi*kT)
-
     def check_velocity_distribution(self, vel, minmax, n_bins, error_tol, kT):
         """check the recorded particle distributions in vel againsta histogram with n_bins bins. Drop velocities outside minmax. Check individual histogram bins up to an accuracy of error_tol agaisnt the analytical result for kT."""
         for i in range(3):
@@ -56,13 +52,13 @@ def check_velocity_distribution(self, vel, minmax, n_bins, error_tol, kT):
             bins = hist[1]
             for j in range(n_bins):
                 found = data[j]
-                expected = self.single_component_maxwell(bins[j], bins[j+1], kT)
+                expected = single_component_maxwell(bins[j], bins[j+1], kT)
                 self.assertLessEqual(abs(found - expected), error_tol)
 
     def test_aa_verify_single_component_maxwell(self):
         """Verifies the normalization of the analytical expression."""
         self.assertLessEqual(
-            abs(self.single_component_maxwell(-10, 10, 4.)-1.), 1E-4)   
+            abs(single_component_maxwell(-10, 10, 4.)-1.), 1E-4)
 
     def check_total_zero(self):
         v_total = np.sum(self.s.part[:].v, axis=0)

testsuite/python/langevin_thermostat.py

@@ -25,6 +25,7 @@
 from time import time
 from espressomd.accumulators import Correlator
 from espressomd.observables import ParticleVelocities, ParticleBodyAngularVelocities
+from tests_common import single_component_maxwell
 
 
 @ut.skipIf(espressomd.has_features("THERMOSTAT_IGNORE_NON_VIRTUAL"),
@@ -45,12 +46,6 @@ class LangevinThermostat(ut.TestCase):
     def setUpClass(cls):
         np.random.seed(42)
 
-    def single_component_maxwell(self, x1, x2, kT):
-        """Integrate the probability density from x1 to x2 using the trapez rule"""
-        x = np.linspace(x1, x2, 1000)
-        return np.trapz(np.exp(-x**2 / (2. * kT)), x) / \
-            np.sqrt(2. * np.pi * kT)
-
     def check_velocity_distribution(self, vel, minmax, n_bins, error_tol, kT):
         """check the recorded particle distributions in vel againsta histogram with n_bins bins. Drop velocities outside minmax. Check individual histogram bins up to an accuracy of error_tol agaisnt the analytical result for kT."""
         for i in range(3):
@@ -60,14 +55,14 @@ def check_velocity_distribution(self, vel, minmax, n_bins, error_tol, kT):
             bins = hist[1]
             for j in range(n_bins):
                 found = data[j]
-                expected = self.single_component_maxwell(
+                expected = single_component_maxwell(
                     bins[j], bins[j + 1], kT)
                 self.assertLessEqual(abs(found - expected), error_tol)
 
     def test_aa_verify_single_component_maxwell(self):
         """Verifies the normalization of the analytical expression."""
         self.assertLessEqual(
-            abs(self.single_component_maxwell(-10, 10, 4.) - 1.), 1E-4)
+            abs(single_component_maxwell(-10, 10, 4.) - 1.), 1E-4)
 
     def test_global_langevin(self):
         """Test for global Langevin parameters."""

testsuite/python/lb_thermostat.py

@@ -0,0 +1,101 @@
+# Copyright (C) 2010-2018 The ESPResSo project
+#
+# This file is part of ESPResSo.
+#
+# ESPResSo is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# ESPResSo is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+import unittest as ut
+import numpy as np
+
+import espressomd.lb
+from tests_common import single_component_maxwell
+
+"""
+Check the Lattice Boltzmann thermostat with respect to the particle velocity distribution.
+
+
+"""
+
+KT = 2.25
+AGRID = 2.5
+VISC = .7
+DENS = 1.7
+TIME_STEP = 0.01
+LB_PARAMS = {'agrid': AGRID,
+             'dens': DENS,
+             'visc': VISC,
+             'fric': 2.0,
+             'tau': TIME_STEP}
+
+
+class LBThermostatCommon(object):
+
+    """Base class of the test that holds the test logic."""
+    lbf = None
+    system = espressomd.System(box_l=[10.0, 10.0, 10.0])
+    system.time_step = TIME_STEP
+    system.cell_system.skin = 0.4 * AGRID
+
+    def prepare(self):
+        self.system.set_random_state_PRNG()
+        self.system.actors.clear()
+        self.system.actors.add(self.lbf)
+        self.system.part.add(
+            pos=np.random.random((250, 3)) * self.system.box_l)
+        self.system.thermostat.set_lb(kT=KT)
+
+    def test_velocity_distribution(self):
+        self.prepare()
+        self.system.integrator.run(100)
+        N = len(self.system.part)
+        loops = 250
+        v_stored = np.zeros((N * loops, 3))
+        for i in range(loops):
+            self.system.integrator.run(10)
+            v_stored[i * N:(i + 1) * N, :] = self.system.part[:].v
+        minmax = 5
+        n_bins = 5
+        error_tol = 0.01
+        for i in range(3):
+            hist = np.histogram(v_stored[:, i], range=(
+                -minmax, minmax), bins=n_bins, normed=False)
+            data = hist[0] / float(v_stored.shape[0])
+            bins = hist[1]
+            for j in range(n_bins):
+                found = data[j]
+                expected = single_component_maxwell(bins[j], bins[j + 1], KT)
+                self.assertLessEqual(abs(found - expected), error_tol)
+
+
+@ut.skipIf(not espressomd.has_features(
+    ['LB']), "Skipping test due to missing features.")
+class LBCPUThermostat(ut.TestCase, LBThermostatCommon):
+
+    """Test for the CPU implementation of the LB."""
+
+    def setUp(self):
+        self.lbf = espressomd.lb.LBFluid(**LB_PARAMS)
+
+
+@ut.skipIf(not espressomd.has_features(
+    ['LB_GPU']), "Skipping test due to missing features.")
+class LBGPUThermostat(ut.TestCase, LBThermostatCommon):
+
+    """Test for the GPU implementation of the LB."""
+
+    def setUp(self):
+        self.lbf = espressomd.lb.LBFluidGPU(**LB_PARAMS)
+
+
+if __name__ == '__main__':
+    ut.main()

testsuite/python/tests_common.py

@@ -588,3 +588,9 @@ class DynamicDict(dict):
     def __getitem__(self, key):
         value = super(DynamicDict, self).__getitem__(key)
         return eval(value, self) if isinstance(value, str) else value
+
+
+def single_component_maxwell(x1, x2, kT):
+    """Integrate the probability density from x1 to x2 using the trapezoidal rule"""
+    x = np.linspace(x1, x2, 1000)
+    return np.trapz(np.exp(-x**2 / (2. * kT)), x) / np.sqrt(2. * np.pi * kT)

testsuite/python/thermalized_bond.py

@@ -23,6 +23,7 @@
 import unittest as ut
 
 import espressomd
+from tests_common import single_component_maxwell
 
 @ut.skipIf(not espressomd.has_features(["MASS"]), "Features not available, skipping test!")
 class ThermalizedBond(ut.TestCase):
@@ -39,13 +40,6 @@ class ThermalizedBond(ut.TestCase):
     def setUpClass(cls):
         np.random.seed(42)
 
-    def single_component_maxwell(self, x1, x2, kT):
-        """Integrate the probability density from x1 to x2 using the trapez rule"""
-
-        x = np.linspace(x1, x2, 1000)
-        return np.trapz(np.exp(-x**2 / (2. * kT)), x) / \
-            np.sqrt(2. * np.pi * kT)
-
     def check_velocity_distribution(self, vel, minmax, n_bins, error_tol, kT):
         """check the recorded particle distributions in vel againsta histogram with n_bins bins. Drop velocities outside minmax. Check individual histogram bins up to an accuracy of error_tol agaisnt the analytical result for kT."""
 
@@ -57,7 +51,7 @@ def check_velocity_distribution(self, vel, minmax, n_bins, error_tol, kT):
 
             for j in range(n_bins):
                 found = data[j]
-                expected = self.single_component_maxwell(
+                expected = single_component_maxwell(
                     bins[j], bins[j + 1], kT)
                 self.assertLessEqual(abs(found - expected), error_tol)