diff --git a/src/script_interface/lbboundaries/LBBoundary.hpp b/src/script_interface/lbboundaries/LBBoundary.hpp
index cfa46fc75a5..55bd59b6198 100644
--- a/src/script_interface/lbboundaries/LBBoundary.hpp
+++ b/src/script_interface/lbboundaries/LBBoundary.hpp
@@ -67,11 +67,9 @@ class LBBoundary : public AutoParameters<LBBoundary> {
     if (name == "get_force") {
       // The get force method uses mpi callbacks on lb cpu
       if (this_node == 0) {
-        const auto rho = lb_lbfluid_get_density();
         const auto agrid = lb_lbfluid_get_agrid();
         const auto tau = lb_lbfluid_get_tau();
-        const double unit_conversion =
-            agrid * agrid * agrid * agrid / rho / tau / tau;
+        const double unit_conversion = agrid / tau / tau;
         return m_lbboundary->get_force() * unit_conversion;
       }
       return none;
diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt
index a06b5e8bf59..1c02b277167 100644
--- a/testsuite/python/CMakeLists.txt
+++ b/testsuite/python/CMakeLists.txt
@@ -121,7 +121,7 @@ python_test(FILE reaction_ensemble.py MAX_NUM_PROC 4)
 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 2 LABELS gpu long)
+python_test(FILE lb_stokes_sphere.py MAX_NUM_PROC 4 LABELS gpu long)
 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)
@@ -179,6 +179,7 @@ python_test(FILE thermalized_bond.py MAX_NUM_PROC 4)
 python_test(FILE thole.py MAX_NUM_PROC 4)
 python_test(FILE lb_switch.py MAX_NUM_PROC 1 LABELS gpu)
 python_test(FILE lb_boundary_velocity.py MAX_NUM_PROC 1)
+python_test(FILE lb_boundary_volume_force.py MAX_NUM_PROC 4)
 python_test(FILE lb_thermo_virtual.py MAX_NUM_PROC 2 LABELS gpu)
 python_test(FILE lb_poiseuille.py MAX_NUM_PROC 4 LABELS gpu)
 python_test(FILE lb_poiseuille_cylinder.py MAX_NUM_PROC 2 LABELS gpu)
diff --git a/testsuite/python/lb_boundary_volume_force.py b/testsuite/python/lb_boundary_volume_force.py
new file mode 100644
index 00000000000..d23ca756d79
--- /dev/null
+++ b/testsuite/python/lb_boundary_volume_force.py
@@ -0,0 +1,116 @@
+# Copyright (C) 2010-2019 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 unittest_decorators as utx
+import numpy as np
+
+import espressomd.lb
+import espressomd.lbboundaries
+import espressomd.shapes
+
+
+"""
+Checks force on lb boundaries for a fluid with a uniform volume force
+"""
+
+
+AGRID = 0.5 
+EXT_FORCE = np.array([-.01, 0.02, 0.03])
+VISC = 3.5
+DENS = 1.5 
+TIME_STEP = 0.05
+LB_PARAMS = {'agrid': AGRID,
+             'dens': DENS,
+             'visc': VISC,
+             'tau': TIME_STEP,
+             'ext_force_density': EXT_FORCE}
+
+
+class LBBoundaryForceCommon:
+
+    """Base class of the test that holds the test logic."""
+    lbf = None
+    system = espressomd.System(box_l=np.array([12.0, 4.0, 4.0]) * AGRID)
+    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
+        accuracy. Then compare the foce balance between force exerted on fluid
+        and forces acting on the boundaries.
+
+        """
+        self.system.actors.clear()
+        self.system.lbboundaries.clear()
+        self.system.actors.add(self.lbf)
+        wall_shape1 = espressomd.shapes.Wall(normal=[1, 0, 0], dist=AGRID)
+        wall_shape2 = espressomd.shapes.Wall(
+            normal=[-1, 0, 0], dist=-(self.system.box_l[0] - AGRID))
+        wall1 = espressomd.lbboundaries.LBBoundary(shape=wall_shape1)
+        wall2 = espressomd.lbboundaries.LBBoundary(shape=wall_shape2)
+
+        self.system.lbboundaries.add(wall1)
+        self.system.lbboundaries.add(wall2)
+        fluid_nodes = self.count_fluid_nodes()
+
+        self.system.integrator.run(20)
+        diff = float("inf")
+        old_val = float("inf")
+        while diff > 0.002:
+            self.system.integrator.run(10)
+            new_val = wall1.get_force()[0]
+            diff = abs(new_val - old_val)
+            old_val = new_val
+
+        surface_area = self.system.box_l[1] * self.system.box_l[2]
+        expected_force = fluid_nodes * AGRID**3 * \
+            np.copy(self.lbf.ext_force_density)
+        measured_force = np.array(wall1.get_force()) + \
+            np.array(wall2.get_force())
+        np.testing.assert_allclose(measured_force, expected_force, atol=2E-2)
+
+
+@utx.skipIfMissingFeatures(['LB_BOUNDARIES', 'EXTERNAL_FORCES'])
+class LBCPUBoundaryForce(ut.TestCase, LBBoundaryForceCommon):
+
+    """Test for the CPU implementation of the LB."""
+
+    def setUp(self):
+        self.lbf = espressomd.lb.LBFluid(**LB_PARAMS)
+
+
+@utx.skipIfMissingGPU()
+@utx.skipIfMissingFeatures(['LB_BOUNDARIES_GPU', 'EXTERNAL_FORCES'])
+class LBGPUBoundaryForce(ut.TestCase, LBBoundaryForceCommon):
+
+    """Test for the GPU implementation of the LB."""
+
+    def setUp(self):
+        self.lbf = espressomd.lb.LBFluidGPU(**LB_PARAMS)
+
+
+if __name__ == '__main__':
+    ut.main()
diff --git a/testsuite/python/lb_shear.py b/testsuite/python/lb_shear.py
index 2673252a5ec..9bc866768c8 100644
--- a/testsuite/python/lb_shear.py
+++ b/testsuite/python/lb_shear.py
@@ -83,7 +83,6 @@ def shear_flow(x, t, nu, v, h, k_max):
 class LBShearCommon:
 
     """Base class of the test that holds the test logic."""
-    lbf = None
     system = espressomd.System(box_l=[H + 2. * AGRID,
                                       W,
                                       W])
@@ -98,9 +97,14 @@ def check_profile(self, shear_plane_normal, shear_direction):
         """
         self.system.lbboundaries.clear()
         self.system.actors.clear()
+        try: 
+            del self.lbf
+        except BaseException:
+            pass
         self.system.box_l = np.max(
             ((W, W, W), shear_plane_normal * (H + 2 * AGRID)), 0)
 
+        self.lbf = self.lb_class(**LB_PARAMS)
         self.system.actors.add(self.lbf)
 
         wall_shape1 = espressomd.shapes.Wall(
@@ -164,8 +168,8 @@ def check_profile(self, shear_plane_normal, shear_direction):
             np.copy(wall1.get_force()),
             -np.copy(wall2.get_force()),
             atol=1E-4)
-        np.testing.assert_allclose(np.copy(wall1.get_force()),
-                                   shear_direction * SHEAR_VELOCITY / H * W**2 * VISC, atol=2E-4)
+        np.testing.assert_allclose(np.dot(np.copy(wall1.get_force()), shear_direction),
+                                   SHEAR_VELOCITY / H * W**2 * dynamic_viscosity, atol=2E-4)
 
     def test(self):
         x = np.array((1, 0, 0), dtype=float)
@@ -185,7 +189,7 @@ class LBCPUShear(ut.TestCase, LBShearCommon):
     """Test for the CPU implementation of the LB."""
 
     def setUp(self):
-        self.lbf = espressomd.lb.LBFluid(**LB_PARAMS)
+        self.lb_class = espressomd.lb.LBFluid
 
 
 @utx.skipIfMissingGPU()
@@ -195,7 +199,7 @@ class LBGPUShear(ut.TestCase, LBShearCommon):
     """Test for the GPU implementation of the LB."""
 
     def setUp(self):
-        self.lbf = espressomd.lb.LBFluidGPU(**LB_PARAMS)
+        self.lb_class = espressomd.lb.LBFluidGPU
 
 
 if __name__ == '__main__':
diff --git a/testsuite/python/lb_stokes_sphere.py b/testsuite/python/lb_stokes_sphere.py
index feb5b90c4e9..103a1d48b0f 100644
--- a/testsuite/python/lb_stokes_sphere.py
+++ b/testsuite/python/lb_stokes_sphere.py
@@ -42,10 +42,10 @@
              'visc': KVISC,
              'tau': TIME_STEP}
 # System setup
-radius = 8 * AGRID
-box_width = 62 * AGRID
+radius = 7 * AGRID
+box_width = 52 * AGRID
 real_width = box_width + 2 * AGRID
-box_length = 62 * AGRID
+box_length = 54 * AGRID
 c_s = np.sqrt(1. / 3. * AGRID**2 / TIME_STEP**2)
 v = [0, 0, 0.2 * c_s]  # The boundary slip
 
@@ -91,10 +91,11 @@ def size(vector):
             return np.sqrt(tmp)
 
         last_force = -1000.
-        stokes_force = 6 * np.pi * KVISC * radius * size(v)
+        dynamic_viscosity = self.lbf.viscosity * self.lbf.density
+        stokes_force = 6 * np.pi * dynamic_viscosity * radius * size(v)
         self.system.integrator.run(35)
         while True:
-            self.system.integrator.run(10)
+            self.system.integrator.run(5)
             force = np.linalg.norm(sphere.get_force())
             if np.abs(last_force - force) < 0.01 * stokes_force:
                 break