Improve testing and code coverage (#4631)

Description of changes:
- improve unit testing of core functionality: P3M, MMM1D, OIF, virtual sites, script interface factory
- remove unreachable code and factor out code duplication
- convert exceptions that are unreachable from the public interface by assertions
- improve Sphinx and Doxygen documentation

.gitlab-ci.yml

@@ -305,6 +305,7 @@ tutorials-samples-maxset:
      with_coverage: 'false'
      with_coverage_python: 'true'
      with_scafacos: 'true'
+     with_stokesian_dynamics: 'true'
      make_check_unit_tests: 'false'
      make_check_python: 'false'
      make_check_tutorials: 'true'

doc/sphinx/analysis.rst

@@ -44,11 +44,11 @@ The different energetic contributions to the total energy can also be obtained (
 
 For example, ::
 
->>> energy = system.analysis.energy()
->>> print(energy["total"])
->>> print(energy["kinetic"])
->>> print(energy["bonded"])
->>> print(energy["non_bonded"])
+    >>> energy = system.analysis.energy()
+    >>> print(energy["total"])
+    >>> print(energy["kinetic"])
+    >>> print(energy["bonded"])
+    >>> print(energy["non_bonded"])
 
 
 .. _Momentum of the system:

doc/sphinx/installation.rst

@@ -721,6 +721,7 @@ The following options control features from external libraries:
 * ``WITH_SCAFACOS``: Build with ScaFaCoS support.
 * ``WITH_GSL``: Build with GSL support.
 * ``WITH_STOKESIAN_DYNAMICS`` Build with Stokesian Dynamics support.
+* ``WITH_PYTHON`` Build with Stokesian Dynamics support.
 
 The following options control code instrumentation:
 
@@ -738,6 +739,7 @@ The following options control how the project is built and tested:
 * ``WITH_CPPCHECK``: Run Cppcheck during compilation.
 * ``WITH_CCACHE``: Enable compiler cache for faster rebuilds.
 * ``WITH_TESTS``: Enable C++ and Python tests.
+* ``WITH_BENCHMARKS``: Enable benchmarks.
 * ``WITH_CUDA_COMPILER`` (string): Select the CUDA compiler.
 * ``CTEST_ARGS`` (string): Arguments passed to the ``ctest`` command.
 * ``TEST_TIMEOUT``: Test timeout.

doc/sphinx/system_setup.rst

@@ -359,7 +359,7 @@ this node is twice as high. For 3 processors, the interactions are 0-0,
 Therefore it is highly recommended that you use N-squared only with an
 odd number of nodes, if with multiple processors at all.
 
-.. _Hybrid:
+.. _Hybrid decomposition:
 
 Hybrid decomposition
 ^^^^^^^^^^^^^^^^^^^^

src/core/MpiCallbacks.hpp

@@ -508,7 +508,7 @@ class MpiCallbacks {
    * @brief Remove callback.
    *
    * Remove the callback id from the callback list.
-   * This is a collective call that must be run on all node.
+   * This is a collective call that must be run on all nodes.
    *
    * @param id Identifier of the callback to remove.
    */
@@ -539,10 +539,8 @@ class MpiCallbacks {
       throw std::logic_error("Callbacks can only be invoked on rank 0.");
     }
 
-    /* Check if callback exists */
-    if (m_callback_map.find(id) == m_callback_map.end()) {
-      throw std::out_of_range("Callback does not exist.");
-    }
+    assert(m_callback_map.find(id) != m_callback_map.end() &&
+           "m_callback_map and m_func_ptr_to_id disagree");
 
     /* Send request to worker nodes */
     boost::mpi::packed_oarchive oa(m_comm);

src/core/electrostatics/specfunc.cpp

@@ -217,30 +217,30 @@ double hzeta(double s, double q) {
   constexpr auto jmax = 12;
   constexpr auto kmax = 10;
 
-  if ((s > max_bits && q < 1.0) || (s > 0.5 * max_bits && q < 0.25))
-    return pow(q, -s);
-  if (s > 0.5 * max_bits && q < 1.0) {
-    double p1 = pow(q, -s);
-    double p2 = pow(q / (1.0 + q), s);
-    double p3 = pow(q / (2.0 + q), s);
+  if ((s > max_bits and q < 1.0) or (s > 0.5 * max_bits and q < 0.25))
+    return std::pow(q, -s);
+  if (s > 0.5 * max_bits and q < 1.0) {
+    auto const p1 = std::pow(q, -s);
+    auto const p2 = std::pow(q / (1.0 + q), s);
+    auto const p3 = std::pow(q / (2.0 + q), s);
     return p1 * (1.0 + p2 + p3);
   }
   /** Euler-Maclaurin summation formula from @cite moshier89a p. 400, with
    *  several typo corrections.
    */
   auto const kmax_q = static_cast<double>(kmax) + q;
-  auto const pmax = pow(kmax_q, -s);
+  auto const pmax = std::pow(kmax_q, -s);
   auto scp = s;
   auto pcp = pmax / kmax_q;
   auto ans = pmax * (kmax_q / (s - 1.0) + 0.5);
 
   for (int k = 0; k < kmax; k++)
-    ans += pow(static_cast<double>(k) + q, -s);
+    ans += std::pow(static_cast<double>(k) + q, -s);
 
   for (int j = 0; j <= jmax; j++) {
     auto const delta = hzeta_c[j + 1] * scp * pcp;
     ans += delta;
-    scp *= (s + 2 * j + 1) * (s + 2 * j + 2);
+    scp *= (s + 2. * j + 1.) * (s + 2. * j + 2.);
     pcp /= Utils::sqr(static_cast<double>(kmax) + q);
   }
 
@@ -251,24 +251,24 @@ double K0(double x) {
   if (x <= 2.0) {
     auto const c = evaluateAsChebychevSeriesAt(bk0_cs, 0.5 * x * x - 1.0);
     auto const i0 = evaluateAsChebychevSeriesAt(bi0_cs, x * x / 4.5 - 1.0);
-    return (-log(x) + Utils::ln_2()) * i0 + c;
+    return (-std::log(x) + Utils::ln_2()) * i0 + c;
   }
   auto const c =
       (x <= 8.0) ? evaluateAsChebychevSeriesAt(ak0_cs, (16.0 / x - 5.0) / 3.0)
                  : evaluateAsChebychevSeriesAt(ak02_cs, 16.0 / x - 1.0);
-  return exp(-x) * c / sqrt(x);
+  return std::exp(-x) * c / std::sqrt(x);
 }
 
 double K1(double x) {
   if (x <= 2.0) {
     auto const c = evaluateAsChebychevSeriesAt(bk1_cs, 0.5 * x * x - 1.0);
     auto const i1 = x * evaluateAsChebychevSeriesAt(bi1_cs, x * x / 4.5 - 1.0);
-    return (log(x) - Utils::ln_2()) * i1 + c / x;
+    return (std::log(x) - Utils::ln_2()) * i1 + c / x;
   }
   auto const c =
       (x <= 8.0) ? evaluateAsChebychevSeriesAt(ak1_cs, (16.0 / x - 5.0) / 3.0)
                  : evaluateAsChebychevSeriesAt(ak12_cs, 16.0 / x - 1.0);
-  return exp(-x) * c / sqrt(x);
+  return std::exp(-x) * c / std::sqrt(x);
 }
 
 /***********************************************************
@@ -287,18 +287,19 @@ static int ak01_orders[] = {
 
 double LPK0(double x) {
   if (x >= 27.) {
-    auto const tmp = .5 * exp(-x) / sqrt(x);
+    auto const tmp = .5 * std::exp(-x) / std::sqrt(x);
     return tmp * ak0_cs[0];
   }
   if (x >= 23.) {
-    auto const tmp = exp(-x) / sqrt(x), xx = (16. / 3.) / x - 5. / 3.;
+    auto const tmp = std::exp(-x) / std::sqrt(x);
+    auto const xx = (16. / 3.) / x - 5. / 3.;
     return tmp * (xx * ak0_cs[1] + 0.5 * ak0_cs[0]);
   }
-  if (x > 2) {
+  if (x > 2.) {
     int j = ak01_orders[static_cast<int>(x) - 2];
     double x2;
     double *s0;
-    if (x <= 8) {
+    if (x <= 8.) {
       s0 = ak0_cs;
       x2 = (2. * 16. / 3.) / x - 2. * 5. / 3.;
     } else {
@@ -312,12 +313,12 @@ double LPK0(double x) {
       d0 = x2 * d0 - dd0 + s0[j];
       dd0 = tmp0;
     }
-    auto const tmp = exp(-x) / sqrt(x);
+    auto const tmp = std::exp(-x) / std::sqrt(x);
     return tmp * (0.5 * (s0[0] + x2 * d0) - dd0);
   }
   /* x <= 2 */
   {
-    /* I0/1 series */
+    /* I0/I1 series */
     int j = 10;
     auto x2 = (2. / 4.5) * x * x - 2.;
     auto dd0 = bi0_cs[j];
@@ -327,7 +328,7 @@ double LPK0(double x) {
       d0 = x2 * d0 - dd0 + bi0_cs[j];
       dd0 = tmp0;
     }
-    auto const tmp = log(x) - Utils::ln_2();
+    auto const tmp = std::log(x) - Utils::ln_2();
     auto const ret = -tmp * (0.5 * (bi0_cs[0] + x2 * d0) - dd0);
 
     /* K0/K1 correction */
@@ -346,11 +347,12 @@ double LPK0(double x) {
 
 double LPK1(double x) {
   if (x >= 27.) {
-    auto const tmp = .5 * exp(-x) / sqrt(x);
+    auto const tmp = .5 * std::exp(-x) / std::sqrt(x);
     return tmp * ak1_cs[0];
   }
   if (x >= 23.) {
-    auto const tmp = exp(-x) / sqrt(x), xx = (16. / 3.) / x - 5. / 3.;
+    auto const tmp = std::exp(-x) / std::sqrt(x);
+    auto const xx = (16. / 3.) / x - 5. / 3.;
     return tmp * (xx * ak1_cs[1] + 0.5 * ak1_cs[0]);
   }
   if (x > 2) {
@@ -371,12 +373,12 @@ double LPK1(double x) {
       d1 = x2 * d1 - dd1 + s1[j];
       dd1 = tmp1;
     }
-    auto const tmp = exp(-x) / sqrt(x);
+    auto const tmp = std::exp(-x) / std::sqrt(x);
     return tmp * (0.5 * (s1[0] + x2 * d1) - dd1);
   }
   /* x <= 2 */
   {
-    /* I0/1 series */
+    /* I0/I1 series */
     int j = 10;
     auto x2 = (2. / 4.5) * x * x - 2.;
     auto dd1 = bi1_cs[j];
@@ -386,7 +388,7 @@ double LPK1(double x) {
       d1 = x2 * d1 - dd1 + bi1_cs[j];
       dd1 = tmp1;
     }
-    auto const tmp = log(x) - Utils::ln_2();
+    auto const tmp = std::log(x) - Utils::ln_2();
     auto const ret = x * tmp * (0.5 * (bi1_cs[0] + x2 * d1) - dd1);
 
     /* K0/K1 correction */
@@ -405,13 +407,14 @@ double LPK1(double x) {
 
 std::pair<double, double> LPK01(double x) {
   if (x >= 27.) {
-    auto const tmp = .5 * exp(-x) / sqrt(x);
+    auto const tmp = .5 * std::exp(-x) / std::sqrt(x);
     auto const k0 = tmp * ak0_cs[0];
     auto const k1 = tmp * ak1_cs[0];
     return {k0, k1};
   }
   if (x >= 23.) {
-    auto const tmp = exp(-x) / sqrt(x), xx = (16. / 3.) / x - 5. / 3.;
+    auto const tmp = std::exp(-x) / std::sqrt(x);
+    auto const xx = (16. / 3.) / x - 5. / 3.;
     auto const k0 = tmp * (xx * ak0_cs[1] + 0.5 * ak0_cs[0]);
     auto const k1 = tmp * (xx * ak1_cs[1] + 0.5 * ak1_cs[0]);
     return {k0, k1};
@@ -440,14 +443,14 @@ std::pair<double, double> LPK01(double x) {
       dd0 = tmp0;
       dd1 = tmp1;
     }
-    auto const tmp = exp(-x) / sqrt(x);
+    auto const tmp = std::exp(-x) / std::sqrt(x);
     auto const k0 = tmp * (0.5 * (s0[0] + x2 * d0) - dd0);
     auto const k1 = tmp * (0.5 * (s1[0] + x2 * d1) - dd1);
     return {k0, k1};
   }
   /* x <= 2 */
   {
-    /* I0/1 series */
+    /* I0/I1 series */
     int j = 10;
     auto x2 = (2. / 4.5) * x * x - 2.;
     auto dd0 = bi0_cs[j];
@@ -461,7 +464,7 @@ std::pair<double, double> LPK01(double x) {
       dd0 = tmp0;
       dd1 = tmp1;
     }
-    auto const tmp = log(x) - Utils::ln_2();
+    auto const tmp = std::log(x) - Utils::ln_2();
     auto k0 = -tmp * (0.5 * (bi0_cs[0] + x2 * d0) - dd0);
     auto k1 = x * tmp * (0.5 * (bi1_cs[0] + x2 * d1) - dd1);
 

src/core/electrostatics/specfunc.hpp

@@ -21,8 +21,8 @@
 
 /** @file
  *  This file contains implementations for some special functions which are
- *  needed by the MMM family of algorithms. This are the modified Hurwitz zeta
- *  function and the modified Bessel functions of first and second kind. The
+ *  needed by the MMM family of algorithms. These are the modified Hurwitz
+ *  zeta function and the modified Bessel functions of second kind. The
  *  implementations are based on the GSL code (see @ref specfunc.cpp for the
  *  original GSL header).
  *
@@ -47,6 +47,7 @@ double hzeta(double order, double x);
 
 /** Modified Bessel function of second kind, order 0. This function was taken
  *  from the GSL code. Precise roughly up to machine precision.
+ *  It is 16 times faster than <tt>std::cyl_bessel_k</tt>.
  *  If @c MMM1D_MACHINE_PREC is not defined, @ref LPK0 is used instead.
  */
 double K0(double x);
@@ -57,21 +58,29 @@ double K0(double x);
  */
 double K1(double x);
 
-/** Bessel function of second kind, order 0, low precision.
+/** Modified Bessel function of second kind, order 0, low precision.
  *  The implementation has an absolute precision of around 10^(-14), which is
  *  comparable to the relative precision sqrt implementation of current
- *  hardware.
+ *  hardware in the ranges @f$ ]0, 8[ @f$ and @f$ ]8, 23[ @f$. Above 23,
+ *  the precision starts to degrade, and above 27 the result drifts and
+ *  slowly converges to 96% of the real value.
+ *  It is 25 times faster than <tt>std::cyl_bessel_k</tt>
+ *  and 1.5 times faster than @ref K0.
  */
 double LPK0(double x);
 
-/** Bessel function of second kind, order 1, low precision.
+/** Modified Bessel function of second kind, order 1, low precision.
  *  The implementation has an absolute precision of around 10^(-14), which is
  *  comparable to the relative precision sqrt implementation of current
- *  hardware.
+ *  hardware in the ranges @f$ ]0, 8[ @f$ and @f$ ]8, 23[ @f$. Above 23,
+ *  the precision starts to degrade, and above 27 the result drifts and
+ *  slowly converges to 111% of the real value.
+ *  It is 25 times faster than <tt>std::cyl_bessel_k</tt>
+ *  and 1.5 times faster than @ref K1.
  */
 double LPK1(double x);
 
-/** Bessel functions of second kind, order 0 and order 1, low precision.
+/** Modified Bessel functions of second kind, order 0 and 1, low precision.
  *  The implementation has an absolute precision of around 10^(-14), which is
  *  comparable to the relative precision sqrt implementation of current
  *  hardware.
@@ -85,8 +94,8 @@ inline double evaluateAsTaylorSeriesAt(Utils::Span<const double> series,
                                        double x) {
   assert(not series.empty());
   auto cnt = static_cast<int>(series.size()) - 1;
-  const double *c = series.data();
-  double r = c[cnt];
+  auto const *c = series.data();
+  auto r = c[cnt];
   while (--cnt >= 0)
     r = r * x + c[cnt];
   return r;
@@ -99,10 +108,10 @@ inline double evaluateAsChebychevSeriesAt(Utils::Span<const double> series,
                                           double x) {
   assert(series.size() >= 3);
 
-  const double *c = series.data();
-  double const x2 = 2.0 * x;
-  double dd = c[series.size() - 1];
-  double d = x2 * dd + c[series.size() - 2];
+  auto const *c = series.data();
+  auto const x2 = 2.0 * x;
+  auto dd = c[series.size() - 1];
+  auto d = x2 * dd + c[series.size() - 2];
   for (auto j = static_cast<int>(series.size()) - 3; j >= 1; j--) {
     auto const tmp = d;
     d = x2 * d - dd + c[j];