Additional LB test cases

testsuite/python/CMakeLists.txt

@@ -192,6 +192,8 @@ python_test(FILE lb_boundary.py MAX_NUM_PROC 2 LABELS gpu)
 python_test(FILE lb_streaming.py MAX_NUM_PROC 4 LABELS gpu)
 python_test(FILE lb_shear.py MAX_NUM_PROC 2 LABELS gpu)
 python_test(FILE lb_thermostat.py MAX_NUM_PROC 2 LABELS gpu)
+python_test(FILE lb_buoyancy_force.py MAX_NUM_PROC 4 LABELS gpu)
+python_test(FILE lb_momentum_conservation.py MAX_NUM_PROC 4 LABELS gpu)
 python_test(FILE p3m_electrostatic_pressure.py MAX_NUM_PROC 2)
 python_test(FILE sigint.py DEPENDENCIES sigint_child.py MAX_NUM_PROC 1)
 python_test(FILE lb_density.py MAX_NUM_PROC 1)

testsuite/python/lb_boundary_volume_force.py

@@ -22,6 +22,8 @@
 import espressomd.lbboundaries
 import espressomd.shapes
 
+from tests_common import count_fluid_nodes
+
 
 """
 Checks force on lb boundaries for a fluid with a uniform volume force
@@ -48,14 +50,6 @@ class LBBoundaryForceCommon:
     system.time_step = TIME_STEP
     system.cell_system.skin = 0.4 * AGRID
 
-    def count_fluid_nodes(self):
-        # Count non-boundary nodes:
-        fluid_nodes = 0
-        for n in self.lbf.nodes():
-            if not n.boundary: 
-                fluid_nodes += 1
-        return fluid_nodes
-
     def test(self):
         """
         Integrate the LB fluid until steady state is reached within a certain
@@ -74,7 +68,7 @@ def test(self):
 
         self.system.lbboundaries.add(wall1)
         self.system.lbboundaries.add(wall2)
-        fluid_nodes = self.count_fluid_nodes()
+        fluid_nodes = count_fluid_nodes(self.lbf)
 
         self.system.integrator.run(20)
         diff = float("inf")

testsuite/python/lb_buoyancy_force.py

@@ -0,0 +1,128 @@
+# 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 espressomd
+from espressomd import lb, lbboundaries, shapes, has_features
+import unittest as ut
+import numpy as np
+import unittest_decorators as utx
+
+from tests_common import count_fluid_nodes
+
+# Define the LB Parameters
+TIME_STEP = 0.01
+AGRID = 0.5 
+KVISC = 6 
+DENS = 2
+G = 0.08
+BOX_SIZE = 24 * AGRID
+
+LB_PARAMS = {'agrid': AGRID,
+             'dens': DENS,
+             'visc': KVISC,
+             'tau': TIME_STEP,
+             'ext_force_density': [0, DENS * G, 0]}
+# System setup
+RADIUS = 8 * AGRID 
+
+
+class Buoyancy(object):
+    """
+    Tests buoyancy force on a sphere in a closed box of lb fluid and 
+    the overall force balance
+
+    """
+    lbf = None
+    system = espressomd.System(box_l=[BOX_SIZE] * 3)
+    system.time_step = TIME_STEP
+    system.cell_system.skin = 0.01
+
+    def test(self):
+        self.system.actors.clear()
+        self.system.lbboundaries.clear()
+        self.system.actors.add(self.lbf)
+
+        # Setup walls
+        for i in range(3):
+            n = np.zeros(3)
+            n[i] = 1
+            self.system.lbboundaries.add(
+                lbboundaries.LBBoundary(shape=shapes.Wall(
+                                        normal=-n, dist=-(self.system.box_l[i] - AGRID))))
+
+            self.system.lbboundaries.add(lbboundaries.LBBoundary(
+                                         shape=shapes.Wall(
+                                             normal=n, dist=AGRID)))
+
+        # setup sphere without slip in the middle
+        sphere = lbboundaries.LBBoundary(shape=shapes.Sphere(
+            radius=RADIUS, center=self.system.box_l / 2, direction=1))
+
+        self.system.lbboundaries.add(sphere)
+
+        sphere_volume = 4. / 3. * np.pi * RADIUS**3
+
+        # Equilibration
+        last_force = -999999
+        self.system.integrator.run(10)
+        while True:
+            self.system.integrator.run(10)
+            force = np.linalg.norm(sphere.get_force())
+
+            if np.linalg.norm(force - last_force) < 0.01:
+                break
+            last_force = force
+
+        # Check force balance
+        boundary_force = np.zeros(3)
+        for b in self.system.lbboundaries:
+            boundary_force += b.get_force()
+
+        fluid_nodes = count_fluid_nodes(self.lbf) 
+        fluid_volume = fluid_nodes * AGRID**3
+        applied_force = fluid_volume * np.array(LB_PARAMS['ext_force_density'])
+
+        np.testing.assert_allclose(
+            boundary_force,
+            applied_force,
+            atol=0.08 * np.linalg.norm(applied_force))
+
+        # Check buoyancy force on the sphere
+        expected_force = np.array(
+            [0, -sphere_volume * DENS * G, 0])
+        np.testing.assert_allclose(
+            np.copy(sphere.get_force()), expected_force,
+            atol=np.linalg.norm(expected_force) * 0.02)
+
+
+@utx.skipIfMissingGPU()
+@utx.skipIfMissingFeatures(["LB_BOUNDARIES_GPU", "EXTERNAL_FORCES"])
+class LBGPUBuoyancy(ut.TestCase, Buoyancy):
+
+    def setUp(self):
+        self.lbf = espressomd.lb.LBFluidGPU(**LB_PARAMS)
+
+
+@utx.skipIfMissingFeatures(["LB_BOUNDARIES", "EXTERNAL_FORCES"])
+class LBCPUBuoyancy(ut.TestCase, Buoyancy):
+
+    def setUp(self):
+        self.lbf = espressomd.lb.LBFluid(**LB_PARAMS)
+
+
+if __name__ == "__main__":
+    ut.main()

testsuite/python/lb_momentum_conservation.py

@@ -0,0 +1,103 @@
+# 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 espressomd
+from espressomd import lb, lbboundaries, shapes, has_features
+import unittest as ut
+import numpy as np
+import sys
+
+# Define the LB Parameters
+TIME_STEP = 0.1
+AGRID = 1.0 
+KVISC = 5
+DENS = 1 
+BOX_SIZE = 6 * AGRID
+F = 1. / BOX_SIZE**3
+GAMMA = 15
+
+LB_PARAMS = {'agrid': AGRID,
+             'dens': DENS,
+             'visc': KVISC,
+             'tau': TIME_STEP,
+             'ext_force_density': [0, F, 0]}
+
+
+class Momentum(object):
+    """
+    Tests momentum conservatoin for an LB coupled to a particle, where opposing
+    forces are applied to LB and particle. The test should uncover issues
+    with boundary and ghost layer handling.
+
+    """
+    lbf = None
+    system = espressomd.System(box_l=[BOX_SIZE] * 3)
+    system.time_step = TIME_STEP
+    system.cell_system.skin = 0.01
+
+    def test(self):
+        self.system.actors.clear()
+        self.system.part.clear()
+        self.system.actors.add(self.lbf)
+        self.system.thermostat.set_lb(LB_fluid=self.lbf, gamma=GAMMA, seed=1)
+
+        applied_force = self.system.volume() * np.array(
+            LB_PARAMS['ext_force_density'])
+        p = self.system.part.add(
+            pos=(0, 0, 0), ext_force=-applied_force, v=[.1, .2, .3])
+
+        # Reach steady state
+        self.system.integrator.run(500)
+        v_final = np.copy(p.v)
+        momentum = self.system.analysis.linear_momentum()
+
+        for i in range(10):
+            self.system.integrator.run(50)
+            # check that momentum stays constant
+            np.testing.assert_allclose(
+                self.system.analysis.linear_momentum(), momentum, atol=2E-4)
+
+            # Check that particle velocity is stationary
+            # up to the acceleration of 1/2 time step
+            np.testing.assert_allclose(np.copy(p.v), v_final, atol=2.2E-3)
+
+        # Make sure, the particle has crossed the periodic boundaries
+        self.assertGreater(
+            np.amax(
+                np.abs(v_final) *
+                self.system.time),
+            BOX_SIZE)
+
+
+@ut.skipIf(not espressomd.gpu_available() or not espressomd.has_features(
+    ['EXTERNAL_FORCES']), "Skipping test due to missing features.")
+class LBGPUMomentum(ut.TestCase, Momentum):
+
+    def setUp(self):
+        self.lbf = espressomd.lb.LBFluidGPU(**LB_PARAMS)
+
+
+@ut.skipIf(not espressomd.has_features(
+    ['EXTERNAL_FORCES']), "Skipping test due to missing features.")
+class LBCPUMomentum(ut.TestCase, Momentum):
+
+    def setUp(self):
+        self.lbf = espressomd.lb.LBFluid(**LB_PARAMS)
+
+
+if __name__ == "__main__":
+    ut.main()

testsuite/python/tests_common.py

@@ -614,3 +614,14 @@ def single_component_maxwell(x1, x2, kT):
 
 def lists_contain_same_elements(list1, list2):
     return len(list1) == len(list2) and sorted(list1) == sorted(list2)
+
+
+def count_fluid_nodes(lbf):
+    """Counts the non-boundary nodes in the passed lb fluid instance."""
+
+    fluid_nodes = 0
+    for n in lbf.nodes():
+        if not n.boundary: 
+            fluid_nodes += 1
+
+    return fluid_nodes