Merge #3287

3287: LB Boundary force: fix unit conversion and tests r=fweik a=RudolfWeeber

Reference solutions for forces on shear planes and for stokes drag need dynamic viscosity, but LB uses kinetmatic viscosity.
This switches the tests to dynamic viscosity and removes a 1/rho from the unit conversion in the script interface.
Also, a force balance test is added, which does not depend on the density.

Fixes #

Description of changes:
 - 


PR Checklist
------------
 - [ ] Tests?
   - [ ] Interface
   - [ ] Core 
 - [ ] Docs?


Co-authored-by: Rudolf Weeber <weeber@icp.uni-stuttgart.de>

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;

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)

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()

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])
@@ -101,6 +100,7 @@ def check_profile(self, shear_plane_normal, shear_direction):
  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 +164,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 +185,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 +195,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__':

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