Additional LB test cases

testsuite/python/CMakeLists.txt

@@ -189,6 +189,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 2 LABELS gpu)
+python_test(FILE lb_momentum_conservation.py MAX_NUM_PROC 2 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_buoyancy_force.py

@@ -0,0 +1,117 @@
+# 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.01
+AGRID = 0.5 
+KVISC = 6 
+DENS = 2
+G = 0.08
+BOX_SIZE = 27 * 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):
+    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)
+
+        expected_force = np.array(
+            [0, -4. / 3. * np.pi * radius**3 * DENS * G, 0])
+        last_force = -999999
+        self.system.integrator.run(100)
+        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()
+        applied_force = ((BOX_SIZE - AGRID)**3 - 4. / 3. * np.pi * radius**3) * \
+            np.array(LB_PARAMS['ext_force_density'])
+        np.testing.assert_allclose(
+            boundary_force,
+            applied_force,
+            rtol=6E-2,
+            atol=0.08 * np.linalg.norm(applied_force))
+        force = np.copy(sphere.get_force())
+        np.testing.assert_allclose(
+            force, expected_force,
+            rtol=0.03,
+            atol=np.linalg.norm(expected_force) * 0.03)
+
+
+@ut.skipIf(not espressomd.gpu_available() or not espressomd.has_features(
+    ['LB_GPU', 'LB_BOUNDARIES_GPU', 'EXTERNAL_FORCES']), "Skipping test due to missing features.")
+class LBGPUBuoyancy(ut.TestCase, Buoyancy):
+
+    def setUp(self):
+        self.lbf = espressomd.lb.LBFluidGPU(**LB_PARAMS)
+
+
+@ut.skipIf(not espressomd.has_features(
+    ['LB', 'LB_BOUNDARIES', 'EXTERNAL_FORCES']), "Skipping test due to missing features.")
+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,95 @@
+# 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.01
+AGRID = 0.5 
+KVISC = 6 
+DENS = 0.5
+
+F = 0.013
+BOX_SIZE = 10 * AGRID
+
+LB_PARAMS = {'agrid': AGRID,
+             'dens': DENS,
+             'visc': KVISC,
+             'tau': TIME_STEP,
+             'ext_force_density': [0, F, 0]}
+
+
+class Momentum(object):
+    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=1, 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)
+
+        # Reach steady state
+        self.system.integrator.run(700)
+        v_final = np.copy(p.v)
+
+        for i in range(3):
+            self.system.integrator.run(100)
+
+            # Check that particle momentum =-fluid momenum
+            # up to the momentum trnasferred in 1/2 time step
+            np.testing.assert_allclose(
+                p.v * p.mass, -
+                    np.array(
+                        self.system.analysis.analyze_linear_momentum(
+                            include_particles=False)),
+                atol=np.linalg.norm(applied_force) * TIME_STEP * 0.55)
+
+            # 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=np.linalg.norm(applied_force) / p.mass * TIME_STEP * 0.55)
+
+
+@ut.skipIf(not espressomd.gpu_available() or not espressomd.has_features(
+    ['LB_GPU', 'LB_BOUNDARIES_GPU', '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(
+    ['LB', 'LB_BOUNDARIES', '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()