tests: Check IBM_VolCons relaxation constant

Make sure volume conservation relaxes the mesh volume towards the
reference volume. Add a callback to query the current volume.

src/core/immersed_boundary/ImmersedBoundaries.hpp

@@ -41,6 +41,11 @@ class ImmersedBoundaries {
       VolumesCurrent.resize(new_size);
     }
   }
+  double get_current_volume(int softID) const {
+    assert(softID >= 0);
+    assert(softID < VolumesCurrent.size());
+    return VolumesCurrent[softID];
+  }
 
 private:
   void calc_volumes(CellStructure &cs);

src/python/espressomd/interactions.pxd

@@ -528,7 +528,7 @@ cdef extern from "thermostat.hpp":
 
 cdef extern from "immersed_boundary/ImmersedBoundaries.hpp":
     cppclass ImmersedBoundaries:
-        pass
+        double get_current_volume(int softID)
 
 cdef extern from "immersed_boundaries.hpp":
     extern ImmersedBoundaries immersed_boundaries

src/python/espressomd/interactions.pyx

@@ -2988,6 +2988,14 @@ class IBM_VolCons(BondedInteraction):
         set_bonded_ia_params(self._bond_id, CoreIBMVolCons(
             self._params["softID"], self._params["kappaV"]))
 
+    def current_volume(self):
+        """
+        Query the current volume of the soft object associated to this bond.
+        The volume is initialized once all :class:`IBM_Triel` bonds have
+        been added and the forces have been recalculated.
+        """
+        return immersed_boundaries.get_current_volume(self._params['softID'])
+
 
 class OifGlobalForces(BondedInteraction):
 

testsuite/python/ibm.py

@@ -22,6 +22,7 @@
 import unittest_decorators as utx
 
 import espressomd
+import espressomd.lb
 import espressomd.interactions
 try:
     from espressomd.virtual_sites import VirtualSitesInertialessTracers
@@ -33,11 +34,12 @@
 class IBM(ut.TestCase):
     '''Test IBM implementation with a Langevin thermostat.'''
     system = espressomd.System(box_l=3 * [10.])
-    system.time_step = 0.05
+    system.time_step = 0.06
     system.cell_system.skin = 0.1
 
     def tearDown(self):
         self.system.part.clear()
+        self.system.actors.clear()
         self.system.thermostat.turn_off()
 
     def compute_dihedral_angle(self, pos0, pos1, pos2, pos3):
@@ -146,7 +148,8 @@ def test_volcons(self):
         positions = list(itertools.product((0, 1), repeat=3))
         positions = positions[:4] + positions[6:] + positions[4:6]
         positions = np.array(positions) - 0.5
-        system.part.add(pos=positions + system.box_l / 2., virtual=8 * [True])
+        mesh_center_ref = np.copy(system.box_l) / 2.
+        system.part.add(pos=positions + mesh_center_ref, virtual=8 * [True])
 
         # Divide the cube. All triangle normals must point inside the mesh.
         # Use the right hand rule to determine the order of the indices.
@@ -172,13 +175,17 @@ def test_volcons(self):
 
         # Run the integrator to initialize the mesh reference volume.
         system.integrator.run(0, recalc_forces=True)
+        self.assertAlmostEqual(volCons.current_volume(), 1., delta=1e-10)
 
         # The restorative force is zero at the moment.
         np.testing.assert_almost_equal(np.copy(self.system.part[:].f), 0.)
 
         # Double the cube dimensions. The volume increases by a factor of 8.
-        system.part[:].pos = 2. * positions + system.box_l / 2.
+        system.part[:].pos = 2. * positions + mesh_center_ref
         system.integrator.run(0, recalc_forces=True)
+        self.assertAlmostEqual(volCons.current_volume(), 8., delta=1e-10)
+
+        # Reference forces for that particular mesh geometry.
         ref_forces = 1.75 * np.array(
             [(1, 2, 2), (2, 1, -2), (2, -1, 1), (1, -2, -1),
              (-1, -2, 2), (-2, -1, -2), (-2, 1, 1), (-1, 2, -1)])
@@ -191,6 +198,50 @@ def test_volcons(self):
         self.assertAlmostEqual(energy['bonded'], 0., delta=1e-10)
         self.assertAlmostEqual(pressure['bonded'], 0., delta=1e-10)
 
+        # Add unthermalized LB.
+        system.thermostat.turn_off()
+        lbf = espressomd.lb.LBFluid(
+            kT=0.0, agrid=2, dens=1, visc=1.8, tau=system.time_step)
+        system.actors.add(lbf)
+        system.thermostat.set_lb(LB_fluid=lbf, act_on_virtual=False, gamma=1)
+
+        # Check the cube is shrinking. The geometrical center shouldn't move.
+        volume_diff_ref = 4.679e-5
+        # warmup
+        system.integrator.run(80)
+        # sampling
+        previous_volume = volCons.current_volume()
+        for _ in range(10):
+            system.integrator.run(20)
+            current_volume = volCons.current_volume()
+            volume_diff = previous_volume - current_volume
+            self.assertLess(current_volume, previous_volume)
+            self.assertAlmostEqual(volume_diff, volume_diff_ref, places=7)
+            previous_volume = current_volume
+            mesh_center = np.mean(self.system.part[:].pos, axis=0)
+            np.testing.assert_allclose(mesh_center, mesh_center_ref)
+
+        # Halve the cube dimensions. The volume decreases by a factor of 8.
+        system.part[:].pos = 0.5 * positions + mesh_center_ref
+        system.integrator.run(0, recalc_forces=True)
+        self.assertAlmostEqual(volCons.current_volume(), 1. / 8., delta=1e-10)
+
+        # Check the cube is expanding. The geometrical center shouldn't move.
+        volume_diff_ref = 7.2118e-7
+        # warmup
+        system.integrator.run(80)
+        # sampling
+        previous_volume = volCons.current_volume()
+        for _ in range(10):
+            system.integrator.run(20)
+            current_volume = volCons.current_volume()
+            volume_diff = current_volume - previous_volume
+            self.assertGreater(current_volume, previous_volume)
+            self.assertAlmostEqual(volume_diff, volume_diff_ref, places=7)
+            previous_volume = current_volume
+            mesh_center = np.mean(self.system.part[:].pos, axis=0)
+            np.testing.assert_allclose(mesh_center, mesh_center_ref)
+
 
 if __name__ == "__main__":
     ut.main()