Merge #3011

3011: tune_skin: add adjust_max_skin feature. r=RudolfWeeber a=pkreissl

Fixes #2961

Description of changes:
 - Add `adjust_max_skin` parameter to `tune_skin()` method, which adds the option to automatically adjust `max_skin` if the passed value is larger than permitted (0.5 * local box length - short range cutoff)
 - Rename core variables to match interface.

As a test you can run the following script adjusted from @RudolfWeeber s minimal script in issue #2961
```python
import espressomd
s = espressomd.System(box_l=[1,1,1])
s.time_step = 0.01
s.non_bonded_inter[0,0].lennard_jones.set_params(epsilon=1,sigma=0.2,cutoff=0.3,shift="auto")
print("=> Tune with adjustment:")
s.cell_system.tune_skin(min_skin=0.1,max_skin=0.6,tol=0.05, int_steps=3, adjust_max_skin=True)
print(s.cell_system.skin)
print("=> and fail without:")
s.cell_system.tune_skin(min_skin=0.1,max_skin=0.6,tol=0.05, int_steps=3)
``` 
resulting in the following output:
```
$ ./pypresso tune_skin.py
=> Tune with adjustment:
0.17500000000000002
=> and fail without:
There were unhandled errors.
ERROR: number of cells 1 is smaller than minimum 8 (interaction range too large or min_num_cells too large)
```



Co-authored-by: Patrick Kreissl <[email protected]>
Co-authored-by: Jean-Noël Grad <[email protected]>

CMakeLists.txt

@@ -235,7 +235,7 @@ if(WITH_CUDA)
   endif()
 endif(WITH_CUDA)
 
-find_package(PythonInterp 3 REQUIRED)
+find_package(PythonInterp 3.2 REQUIRED)
 
 if(WITH_PYTHON)
   find_package(Cython 0.23 REQUIRED)

maintainer/benchmarks/p3m.py

@@ -151,14 +151,16 @@
 # tuning and equilibration
 system.integrator.run(min(3 * measurement_steps, 1000))
 print("Tune skin: {}".format(system.cell_system.tune_skin(
-    min_skin=0.4, max_skin=1.6, tol=0.05, int_steps=100)))
+    min_skin=0.4, max_skin=1.6, tol=0.05, int_steps=100,
+    adjust_max_skin=True)))
 system.integrator.run(min(3 * measurement_steps, 3000))
 print("Tune p3m")
 p3m = electrostatics.P3M(prefactor=args.prefactor, accuracy=1e-4)
 system.actors.add(p3m)
 system.integrator.run(min(3 * measurement_steps, 3000))
 print("Tune skin: {}".format(system.cell_system.tune_skin(
-    min_skin=1.0, max_skin=1.6, tol=0.05, int_steps=100)))
+    min_skin=1.0, max_skin=1.6, tol=0.05, int_steps=100,
+    adjust_max_skin=True)))
 
 
 if not args.visualizer:

src/core/tuning.cpp

@@ -24,12 +24,14 @@
 #include "communication.hpp"
 #include "errorhandling.hpp"
 #include "global.hpp"
+#include "grid.hpp"
 #include "integrate.hpp"
 #include <limits>
 #include <sys/resource.h>
 #include <sys/time.h>
 #include <utils/statistics/RunningAverage.hpp>
 
+#include <boost/range/algorithm/min_element.hpp>
 int timing_samples = 10;
 
 /**
@@ -94,21 +96,28 @@ static double time_calc(int rds) {
   return 1000. * (tock - tick) / rds;
 }
 
-void tune_skin(double min, double max, double tol, int steps) {
+void tune_skin(double min_skin, double max_skin, double tol, int int_steps,
+               bool adjust_max_skin) {
   skin_set = true;
 
-  double a = min;
-  double b = max;
+  double a = min_skin;
+  double b = max_skin;
   double time_a, time_b;
 
+  auto const min_local_box_l = *boost::min_element(local_geo.length());
+  double const max_permissible_skin = 0.5 * min_local_box_l - max_cut;
+
+  if (adjust_max_skin and max_skin > max_permissible_skin)
+    b = max_permissible_skin;
+
   while (fabs(a - b) > tol) {
     skin = a;
     mpi_bcast_parameter(FIELD_SKIN);
-    time_a = time_calc(steps);
+    time_a = time_calc(int_steps);
 
     skin = b;
     mpi_bcast_parameter(FIELD_SKIN);
-    time_b = time_calc(steps);
+    time_b = time_calc(int_steps);
 
     if (time_a > time_b) {
       a = 0.5 * (a + b);

src/core/tuning.hpp

@@ -49,6 +49,7 @@ double diffTime();
 /** sets the optimal \ref skin between min and max by bisection to tolerance
  *  tol.
  */
-void tune_skin(double min, double max, double tol, int steps);
+void tune_skin(double min_skin, double max_skin, double tol, int int_steps,
+               bool adjust_max_skin);
 
 #endif

src/python/espressomd/cellsystem.pxd

@@ -18,6 +18,8 @@
 #
 
 # @TODO: shouldn't these global definitions be used via global_variables?
+
+from libcpp cimport bool
 from libcpp.vector cimport vector
 from libcpp.pair cimport pair
 
@@ -41,4 +43,4 @@ cdef extern from "layered.hpp":
     int n_layers_ "n_layers"
 
 cdef extern from "tuning.hpp":
-    cdef void c_tune_skin "tune_skin" (double min, double max, double tol, int steps)
+    cdef void c_tune_skin "tune_skin" (double min_skin, double max_skin, double tol, int int_steps, bool adjust_max_skin)

src/python/espressomd/cellsystem.pyx

@@ -280,7 +280,7 @@ cdef class CellSystem(object):
             return skin
 
     def tune_skin(self, min_skin=None, max_skin=None, tol=None,
-                  int_steps=None):
+                  int_steps=None, adjust_max_skin=False):
         """
         Tunes the skin by measuring the integration time and bisecting over the
         given range of skins. The best skin is set in the simulation core.
@@ -295,11 +295,16 @@ cdef class CellSystem(object):
                 Accuracy in skin to tune to.
         'int_steps' : :obj:`int`
                       Integration steps to time.
+        'adjust_max_skin' : :obj:`bool`, optional
+                            If ``True``, the value of ``max_skin`` is reduced
+                            to the maximum permissible skin (in case the passed
+                            value is too large). Set to ``False`` by default.
 
         Returns
         -------
         :attr:`skin`
 
         """
-        c_tune_skin(min_skin, max_skin, tol, int_steps)
+        c_tune_skin(min_skin, max_skin, tol, int_steps, adjust_max_skin)
+        handle_errors("Error during tune_skin")
         return self.skin

src/python/espressomd/interactions.pyx

@@ -1793,7 +1793,7 @@ class NonBondedInteractionHandle(object):
 cdef class NonBondedInteractions(object):
     """
     Access to non-bonded interaction parameters via [i,j], where i,j are particle
-    types. Returns NonBondedInteractionHandle.
+    types. Returns a :class:`NonBondedInteractionHandle` object.
     Also: access to force capping.
 
     """

src/python/espressomd/tuning.pxd

@@ -1,3 +1,4 @@
+from libcpp cimport bool
 
 cdef extern from "tuning.hpp":
-    cdef void tune_skin(double min, double max, double tol, int steps)
+    cdef void tune_skin(double min_skin, double max_skin, double tol, int int_steps, bool adjust_max_skin)

testsuite/python/CMakeLists.txt

@@ -77,6 +77,7 @@ foreach(TEST_COMBINATION lb.cpu-p3m.cpu-lj-lbtherm;lb.gpu-p3m.cpu-lj-lvtherm;ek.
   endforeach(TEST_BINARY)
 endforeach(TEST_COMBINATION)
 python_test(FILE cellsystem.py MAX_NUM_PROC 4)
+python_test(FILE tune_skin.py MAX_NUM_PROC 1)
 python_test(FILE constraint_homogeneous_magnetic_field.py MAX_NUM_PROC 4)
 python_test(FILE constraint_shape_based.py MAX_NUM_PROC 2)
 python_test(FILE coulomb_cloud_wall.py MAX_NUM_PROC 4 LABELS gpu)

testsuite/python/polymer.py

@@ -14,7 +14,6 @@
 #
 # 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 sys
 import unittest as ut
 import numpy as np
 import espressomd
@@ -181,9 +180,7 @@ def test_respect_constraints_wall(self):
             self.assertGreaterEqual(z, 0.5 * self.box_l)
 
         # assert that illegal start position raises error
-        assertRaisesRegex = self.assertRaisesRegexp if sys.version_info < (
-            3, 2) else self.assertRaisesRegex
-        with assertRaisesRegex(Exception, 'Invalid start positions.'):
+        with self.assertRaisesRegex(Exception, 'Invalid start positions.'):
             illegal_start = np.array([[1., 1., 0.2 * self.box_l]])
             positions = polymer.positions(
                 n_polymers=1,

testsuite/python/test_checkpoint.py

@@ -15,7 +15,6 @@
 # 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 sys
 import subprocess
 import unittest as ut
 import unittest_decorators as utx
@@ -53,9 +52,7 @@ def test_LB(self):
         cpt_path = self.checkpoint.checkpoint_dir + "/lb{}.cpt"
         with self.assertRaises(RuntimeError):
             lbf.load_checkpoint(cpt_path.format("-corrupted"), cpt_mode)
-        assertRaisesRegex = self.assertRaisesRegexp if sys.version_info < (
-            3, 2) else self.assertRaisesRegex
-        with assertRaisesRegex(RuntimeError, 'grid dimensions mismatch'):
+        with self.assertRaisesRegex(RuntimeError, 'grid dimensions mismatch'):
             lbf.load_checkpoint(cpt_path.format("-wrong-boxdim"), cpt_mode)
         lbf.load_checkpoint(cpt_path.format(""), cpt_mode)
         precision = 9 if "LB.CPU" in modes else 5

testsuite/python/tune_skin.py

@@ -0,0 +1,52 @@
+# 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 unittest as ut
+import unittest_decorators as utx
+import espressomd
+
+
+@utx.skipIfMissingFeatures("LENNARD_JONES")
+class TuneSkin(ut.TestCase):
+    system = espressomd.System(box_l=[1.35, 2.4, 1.7])
+    system.time_step = 0.01
+
+    def setUp(self):
+        self.system.non_bonded_inter[0, 0].lennard_jones.set_params(
+            epsilon=1,
+            sigma=0.2,
+            cutoff=0.3,
+            shift="auto")
+
+    def test_fails_without_adjustment(self):
+        with self.assertRaisesRegex(Exception, 'Error during tune_skin'):
+            self.system.cell_system.tune_skin(
+                min_skin=0.1,
+                max_skin=0.6,
+                tol=0.05,
+                int_steps=3)
+
+    def test_works_with_adjustment(self):
+        self.system.cell_system.tune_skin(
+            min_skin=0.1,
+            max_skin=0.6,
+            tol=0.05,
+            int_steps=3,
+            adjust_max_skin=True)
+
+if __name__ == "__main__":
+    ut.main()