From a30a627492fab9ce742f7e5d7e04898b3d10ff1f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 8 Dec 2022 23:12:01 +0100 Subject: [PATCH 1/6] Modernize, document and unit test special functions --- src/core/electrostatics/specfunc.cpp | 63 ++++++------- src/core/electrostatics/specfunc.hpp | 35 +++++--- src/core/unit_tests/CMakeLists.txt | 7 +- src/core/unit_tests/specfunc_test.cpp | 122 ++++++++++++++++++++++++++ 4 files changed, 183 insertions(+), 44 deletions(-) create mode 100644 src/core/unit_tests/specfunc_test.cpp diff --git a/src/core/electrostatics/specfunc.cpp b/src/core/electrostatics/specfunc.cpp index 79057995562..698995b8129 100644 --- a/src/core/electrostatics/specfunc.cpp +++ b/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(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(k) + q, -s); + ans += std::pow(static_cast(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(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(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 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 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 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); diff --git a/src/core/electrostatics/specfunc.hpp b/src/core/electrostatics/specfunc.hpp index 55df4af2652..8923d3346dc 100644 --- a/src/core/electrostatics/specfunc.hpp +++ b/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 std::cyl_bessel_k. * 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 std::cyl_bessel_k + * 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 std::cyl_bessel_k + * 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 series, double x) { assert(not series.empty()); auto cnt = static_cast(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 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(series.size()) - 3; j >= 1; j--) { auto const tmp = d; d = x2 * d - dd + c[j]; diff --git a/src/core/unit_tests/CMakeLists.txt b/src/core/unit_tests/CMakeLists.txt index 82465c13b89..9be4e36a851 100644 --- a/src/core/unit_tests/CMakeLists.txt +++ b/src/core/unit_tests/CMakeLists.txt @@ -35,7 +35,7 @@ unit_test(NAME MpiCallbacks_test SRC MpiCallbacks_test.cpp DEPENDS espresso::utils Boost::mpi MPI::MPI_CXX NUM_PROC 2) unit_test(NAME ParticleIterator_test SRC ParticleIterator_test.cpp DEPENDS espresso::utils) -unit_test(NAME p3m_test SRC p3m_test.cpp DEPENDS espresso::utils) +unit_test(NAME p3m_test SRC p3m_test.cpp DEPENDS espresso::utils espresso::core) unit_test(NAME link_cell_test SRC link_cell_test.cpp DEPENDS espresso::utils) unit_test(NAME Particle_test SRC Particle_test.cpp DEPENDS espresso::utils Boost::serialization) @@ -70,3 +70,8 @@ unit_test(NAME bonded_interactions_map_test SRC bonded_interactions_map_test.cpp DEPENDS espresso::core) unit_test(NAME bond_breakage_test SRC bond_breakage_test.cpp DEPENDS espresso::core) +if(NOT CMAKE_CXX_COMPILER_ID STREQUAL "AppleClang") + # AppleClang doesn't implement C++17's mathematical special functions + unit_test(NAME specfunc_test SRC specfunc_test.cpp DEPENDS espresso::utils + espresso::core) +endif() diff --git a/src/core/unit_tests/specfunc_test.cpp b/src/core/unit_tests/specfunc_test.cpp new file mode 100644 index 00000000000..6340312f344 --- /dev/null +++ b/src/core/unit_tests/specfunc_test.cpp @@ -0,0 +1,122 @@ +/* + * Copyright (C) 2022 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 . + */ + +#define BOOST_TEST_MODULE special functions test +#define BOOST_TEST_DYN_LINK +#include + +#include "electrostatics/specfunc.hpp" + +#include +#include + +constexpr auto eps = 100. * std::numeric_limits::epsilon(); + +BOOST_AUTO_TEST_CASE(hurwitz_zeta_function) { + constexpr auto max_bits = 54.0; + // test cases where an exact, closed-form expression exists + auto delta = 0.025; + auto x = delta / 2.; + while (x < 0.25) { + auto order = max_bits / 2. + 0.01; + BOOST_TEST_INFO("with parameter x = " << x); + BOOST_CHECK_CLOSE(hzeta(order, x), std::pow(x, -order), eps); + x += delta; + } + x = delta / 2.; + while (x < 1.) { + auto order = max_bits + 0.01; + BOOST_TEST_INFO("with parameter x = " << x); + BOOST_CHECK_CLOSE(hzeta(order, x), std::pow(x, -order), eps); + x += delta; + } + x = delta / 2.; + while (x < 1.) { + auto order = max_bits / 2. + 0.01; + auto ref = std::pow(x, -order); + ref *= (1. + std::pow(x / (1. + x), order) + std::pow(x / (2. + x), order)); + BOOST_TEST_INFO("with parameter x = " << x); + BOOST_CHECK_CLOSE(hzeta(order, x), ref, eps); + x += delta; + } + // Hurwitz zeta is a generalization of the Riemann zeta + for (auto order = 2; order < 128; ++order) { + BOOST_TEST_INFO("with parameter order = " << order); + BOOST_CHECK_CLOSE(hzeta(order, 1.), std::riemann_zeta(order), 10. * eps); + } +} + +BOOST_AUTO_TEST_CASE(bessel_series_high_precision) { + auto delta = 0.02; + auto x = delta; + while (x < 20.) { + BOOST_TEST_INFO("with parameter x = " << x); + BOOST_CHECK_CLOSE(K0(x), std::cyl_bessel_k(0, x), 10. * eps); + BOOST_CHECK_CLOSE(K1(x), std::cyl_bessel_k(1, x), 10. * eps); + x += delta; + } + delta = 0.5; + while (x < 64.) { + BOOST_TEST_INFO("with parameter x = " << x); + BOOST_CHECK_CLOSE(K0(x), std::cyl_bessel_k(0, x), 10. * eps); + BOOST_CHECK_CLOSE(K1(x), std::cyl_bessel_k(1, x), 10. * eps); + x += delta; + } +} + +BOOST_AUTO_TEST_CASE(bessel_series_low_precision) { + auto const check = [](double x, double tol) { + BOOST_TEST_INFO("with parameter x = " << x); + BOOST_CHECK_CLOSE(LPK0(x), std::cyl_bessel_k(0, x), tol); + BOOST_TEST_INFO("with parameter x = " << x); + BOOST_CHECK_CLOSE(LPK1(x), std::cyl_bessel_k(1, x), tol); + BOOST_CHECK_CLOSE(LPK01(x).first, LPK0(x), eps); + BOOST_CHECK_CLOSE(LPK01(x).second, LPK1(x), eps); + }; + auto delta = 0.02; + auto x = delta; + // important breakpoints: x=2, x=8, x=23, x=27 (different interpolation) + while (x < 2.) { + check(x, 20. * eps); + x += delta; + } + while (x < 8.) { + check(x, 1e5 * eps); + x += delta; + } + delta = 0.04; + while (x < 13.) { + check(x, 1e8 * eps); + x += delta; + } + while (x < 23.) { + check(x, 2e10 * eps); + x += delta; + } + delta = 0.08; + while (x < 27.) { + check(x, 0.6); + x += delta; + } + delta = 0.64; + while (x < 64.) { + check(x, 20.); + x += delta; + } +} From 5c69b9b4f7573bb6b42c30e53eb26f27d04e481d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 9 Dec 2022 20:34:54 +0100 Subject: [PATCH 2/6] Document and unit test virtual sites interface --- doc/sphinx/system_setup.rst | 2 +- src/python/espressomd/virtual_sites.py | 13 +++++++++--- testsuite/python/save_checkpoint.py | 2 ++ testsuite/python/test_checkpoint.py | 2 ++ testsuite/python/virtual_sites_relative.py | 24 ++++++++++++++++------ 5 files changed, 33 insertions(+), 10 deletions(-) diff --git a/doc/sphinx/system_setup.rst b/doc/sphinx/system_setup.rst index 5a3837e3f67..a8ab3c23bf2 100644 --- a/doc/sphinx/system_setup.rst +++ b/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 ^^^^^^^^^^^^^^^^^^^^ diff --git a/src/python/espressomd/virtual_sites.py b/src/python/espressomd/virtual_sites.py index ba047c11b44..ab29f03448c 100644 --- a/src/python/espressomd/virtual_sites.py +++ b/src/python/espressomd/virtual_sites.py @@ -30,8 +30,15 @@ class ActiveVirtualSitesHandle(ScriptInterfaceHelper): Attributes ---------- - implementation : - instance of a virtual sites implementation + have_quaternion : :obj:`bool` + Whether the virtual sites has a quaternion (only relevant for + :class:`~espressomd.virtual_sites.VirtualSitesRelative`). + override_cutoff_check : :obj:`bool` + Whether to disable the sanity check that triggers when attempting + to set up a virtual site too far away from the real particle in a + MPI-parallel simulation with more than 1 core. Disabling this + check is not recommended; it is only relevant in rare situations + when using the :ref:`Hybrid decomposition` cell system. """ _so_name = "VirtualSites::ActiveVirtualSitesHandle" @@ -47,7 +54,7 @@ class VirtualSitesOff(ScriptInterfaceHelper): @script_interface_register class VirtualSitesInertialessTracers(ScriptInterfaceHelper): - """Virtual sites which are advected with an lb fluid without inertia. + """Virtual sites which are advected with an LB fluid without inertia. Forces are on them are transferred to the fluid instantly. """ diff --git a/testsuite/python/save_checkpoint.py b/testsuite/python/save_checkpoint.py index 2bea13d10a9..7bcbee17f98 100644 --- a/testsuite/python/save_checkpoint.py +++ b/testsuite/python/save_checkpoint.py @@ -225,6 +225,8 @@ if espressomd.has_features(['VIRTUAL_SITES', 'VIRTUAL_SITES_RELATIVE']): system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative( have_quaternion=True) + system.virtual_sites.have_quaternion = True + system.virtual_sites.override_cutoff_check = True p2.vs_auto_relate_to(p1) # non-bonded interactions diff --git a/testsuite/python/test_checkpoint.py b/testsuite/python/test_checkpoint.py index 12bf3ae6739..35d5f27c790 100644 --- a/testsuite/python/test_checkpoint.py +++ b/testsuite/python/test_checkpoint.py @@ -419,6 +419,8 @@ def test_virtual_sites(self): self.assertIsInstance( system.virtual_sites, espressomd.virtual_sites.VirtualSitesRelative) + self.assertTrue(system.virtual_sites.have_quaternion) + self.assertTrue(system.virtual_sites.override_cutoff_check) def test_mean_variance_calculator(self): np.testing.assert_array_equal( diff --git a/testsuite/python/virtual_sites_relative.py b/testsuite/python/virtual_sites_relative.py index 61f751ef959..271058ce495 100644 --- a/testsuite/python/virtual_sites_relative.py +++ b/testsuite/python/virtual_sites_relative.py @@ -84,12 +84,22 @@ def test_aa_method_switching(self): self.assertIsInstance( self.system.virtual_sites, espressomd.virtual_sites.VirtualSitesOff) + self.assertFalse(self.system.virtual_sites.have_quaternion) + self.assertFalse(self.system.virtual_sites.override_cutoff_check) + + # Set properties + self.system.virtual_sites.have_quaternion = True + self.system.virtual_sites.override_cutoff_check = True + self.assertTrue(self.system.virtual_sites.have_quaternion) + self.assertTrue(self.system.virtual_sites.override_cutoff_check) # Switch implementation self.system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative() self.assertIsInstance( self.system.virtual_sites, espressomd.virtual_sites.VirtualSitesRelative) + self.assertFalse(self.system.virtual_sites.have_quaternion) + self.assertFalse(self.system.virtual_sites.override_cutoff_check) def test_vs_quat(self): self.system.time_step = 0.01 @@ -233,9 +243,11 @@ def test_pos_vel_forces(self): self.assertAlmostEqual(np.linalg.norm(t_exp - t), 0., delta=1E-6) def run_test_lj(self): - """This fills the system with vs-based dumbells, adds a lj potential, - integrates and verifies forces. This is to make sure that no pairs - get lost or are outdated in the short range loop""" + """ + This fills the system with vs-based dumbbells, adds a LJ potential, + integrates and verifies forces. This is to make sure that no pairs + get lost or are outdated in the short range loop. + """ system = self.system system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative() # Parameters @@ -259,8 +271,8 @@ def run_test_lj(self): system.time_step = 0.01 system.thermostat.turn_off() - # Dumbells consist of 2 virtual lj spheres + central particle w/o interactions - # For n spheres, n/2 dumbells. + # Dumbbells consist of 2 virtual lj spheres + central particle + # w/o interactions. For n spheres, n/2 dumbbells. for i in range(n // 2): # Type=1, i.e., no lj ia for the center of mass particles p3i = system.part.add( @@ -382,4 +394,4 @@ def test_zz_pressure_tensor(self): if __name__ == "__main__": - ut.main() + ut.main(verbosity=2) From 3ea0681103a78b7dadc36425d0b79082223d171c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 13 Dec 2022 18:23:35 +0100 Subject: [PATCH 3/6] Fix regressions in user guide and doxygen --- doc/sphinx/analysis.rst | 10 +++++----- doc/sphinx/installation.rst | 2 ++ src/core/analysis/statistics.hpp | 1 + src/core/magnetostatics/dipolar_direct_sum.cpp | 1 + 4 files changed, 9 insertions(+), 5 deletions(-) diff --git a/doc/sphinx/analysis.rst b/doc/sphinx/analysis.rst index a4628ccd921..ae8a9ba0654 100644 --- a/doc/sphinx/analysis.rst +++ b/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: diff --git a/doc/sphinx/installation.rst b/doc/sphinx/installation.rst index 54c55e5050b..58238893ece 100644 --- a/doc/sphinx/installation.rst +++ b/doc/sphinx/installation.rst @@ -753,6 +753,7 @@ The following options control features from external libraries: * ``ESPRESSO_BUILD_WITH_SCAFACOS``: Build with ScaFaCoS support. * ``ESPRESSO_BUILD_WITH_GSL``: Build with GSL support. * ``ESPRESSO_BUILD_WITH_STOKESIAN_DYNAMICS`` Build with Stokesian Dynamics support. +* ``ESPRESSO_BUILD_WITH_PYTHON``: Build with the Python interface. The following options control code instrumentation: @@ -770,6 +771,7 @@ The following options control how the project is built and tested: * ``ESPRESSO_BUILD_WITH_CPPCHECK``: Run Cppcheck during compilation. * ``ESPRESSO_BUILD_WITH_CCACHE``: Enable compiler cache for faster rebuilds. * ``ESPRESSO_BUILD_TESTS``: Enable C++ and Python tests. +* ``ESPRESSO_BUILD_BENCHMARKS``: Enable benchmarks. * ``ESPRESSO_CUDA_COMPILER`` (string): Select the CUDA compiler. * ``ESPRESSO_CTEST_ARGS`` (string): Arguments passed to the ``ctest`` command. * ``ESPRESSO_TEST_TIMEOUT``: Test timeout. diff --git a/src/core/analysis/statistics.hpp b/src/core/analysis/statistics.hpp index b8d5e70076e..3de856bfd96 100644 --- a/src/core/analysis/statistics.hpp +++ b/src/core/analysis/statistics.hpp @@ -98,6 +98,7 @@ std::vector> structure_factor(PartCfg &partCfg, std::vector const &p_types, int order); /** Calculate the center of mass of a special type of the current configuration. + * @param partCfg particle collection * @param p_type type of the particle */ Utils::Vector3d center_of_mass(PartCfg &partCfg, int p_type); diff --git a/src/core/magnetostatics/dipolar_direct_sum.cpp b/src/core/magnetostatics/dipolar_direct_sum.cpp index b668d33529d..63853e7fe57 100644 --- a/src/core/magnetostatics/dipolar_direct_sum.cpp +++ b/src/core/magnetostatics/dipolar_direct_sum.cpp @@ -166,6 +166,7 @@ struct PosMom { * at all. If false, distances are calulated as Euclidean * distances, and not using minimum image convention. * @param ncut Number of replicas in each direction. + * @param box_l Box dimensions. * @param init Initial value of the sum. * @param f Binary operation mapping distance and moment of the * interaction partner to the value to be summed up for this pair. From 04d9211e44dcb1154dc967f8b1d077f8fa3a267e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 15 Dec 2022 20:40:52 +0100 Subject: [PATCH 4/6] Modernize, optimize and unit test OIF code --- src/core/object-in-fluid/oif_local_forces.hpp | 38 ++++------ src/python/object_in_fluid/oif_classes.py | 20 ++--- testsuite/python/oif_volume_conservation.py | 76 ++++++++++++------- 3 files changed, 75 insertions(+), 59 deletions(-) diff --git a/src/core/object-in-fluid/oif_local_forces.hpp b/src/core/object-in-fluid/oif_local_forces.hpp index 18653cf9477..3a07c67595e 100644 --- a/src/core/object-in-fluid/oif_local_forces.hpp +++ b/src/core/object-in-fluid/oif_local_forces.hpp @@ -98,13 +98,6 @@ struct OifLocalForcesBond { } }; -/** @details see eq. (19) in @cite dupin07a */ -inline double KS(double lambda) { - double res; - res = (pow(lambda, 0.5) + pow(lambda, -2.5)) / (lambda + pow(lambda, -3.)); - return res; -} - /** Compute the OIF local forces. * See @cite dupin07a, @cite jancigova16a. * @param p2 %Particle of triangle 1. @@ -126,24 +119,25 @@ OifLocalForcesBond::calc_forces(Particle const &p2, Particle const &p1, Utils::Vector3d force1{}, force2{}, force3{}, force4{}; - // non-linear stretching - if (ks > TINY_OIF_ELASTICITY_COEFFICIENT) { - auto const dx = fp2 - fp3; - auto const len = dx.norm(); - auto const dr = len - r0; - auto const lambda = 1.0 * len / r0; - auto const fac = -ks * KS(lambda) * dr; // no normalization - auto const f = (fac / len) * dx; - force2 += f; - force3 -= f; - } - - // linear stretching - if (kslin > TINY_OIF_ELASTICITY_COEFFICIENT) { + // surface strain constraint + if (ks > TINY_OIF_ELASTICITY_COEFFICIENT or + kslin > TINY_OIF_ELASTICITY_COEFFICIENT) { auto const dx = fp2 - fp3; auto const len = dx.norm(); auto const dr = len - r0; - auto const fac = -kslin * dr; // no normalization + auto fac = 0.; + // linear stretching + if (kslin > TINY_OIF_ELASTICITY_COEFFICIENT) { + fac -= kslin * dr; // no normalization + } + // non-linear stretching + if (ks > TINY_OIF_ELASTICITY_COEFFICIENT) { + /** For non-linear stretching, see eq. (19) in @cite dupin07a */ + auto const lambda = len / r0; + auto const spring = (std::pow(lambda, 0.5) + std::pow(lambda, -2.5)) / + (lambda + std::pow(lambda, -3.)); + fac -= ks * spring * dr; // no normalization + } auto const f = (fac / len) * dx; force2 += f; force3 -= f; diff --git a/src/python/object_in_fluid/oif_classes.py b/src/python/object_in_fluid/oif_classes.py index 50a287635fd..f471f18d45f 100644 --- a/src/python/object_in_fluid/oif_classes.py +++ b/src/python/object_in_fluid/oif_classes.py @@ -933,16 +933,18 @@ def volume(self): return self.mesh.volume() def diameter(self): - max_distance = 0.0 + """ + Compute the maximal pairwise distance between all mesh points. + Computational complexity is :math:`\\mathcal{O}(N^2)`. + """ n_points = len(self.mesh.points) - for i in range(0, n_points): - for j in range(i + 1, n_points): - p1 = self.mesh.points[i].get_pos() - p2 = self.mesh.points[j].get_pos() - tmp_dist = vec_distance(p1, p2) - if tmp_dist > max_distance: - max_distance = tmp_dist - return max_distance + pos = np.array(list(map(lambda p: p.get_pos(), self.mesh.points))) + max_dist_sq = 0.0 + for i in range(0, n_points - 1): + dist_sq = np.max((np.square(pos[i] - pos[i + 1:])).sum(axis=1)) + if dist_sq > max_dist_sq: + max_dist_sq = dist_sq + return np.sqrt(max_dist_sq) def get_n_nodes(self): return self.mesh.get_n_nodes() diff --git a/testsuite/python/oif_volume_conservation.py b/testsuite/python/oif_volume_conservation.py index e437ad67cab..054787f24f0 100644 --- a/testsuite/python/oif_volume_conservation.py +++ b/testsuite/python/oif_volume_conservation.py @@ -22,6 +22,8 @@ import unittest as ut import unittest_decorators as utx import tests_common +import scipy.optimize +import object_in_fluid as oif @utx.skipIfMissingFeatures("MASS") @@ -30,59 +32,77 @@ class OifVolumeConservation(ut.TestCase): """Loads a soft elastic sphere via object_in_fluid, stretches it and checks restoration of original volume due to elastic forces.""" - def test(self): - import object_in_fluid as oif + system = espressomd.System(box_l=(50., 50., 50.)) + system.time_step = 0.4 + system.cell_system.skin = 0.5 - system = espressomd.System(box_l=(10, 10, 10)) - self.assertEqual(system.max_oif_objects, 0) - system.time_step = 0.4 - system.cell_system.skin = 0.5 + def check_relaxation(self, **kwargs): + self.system.part.clear() + self.system.thermostat.turn_off() + half_box_l = np.copy(self.system.box_l) / 2. # creating the template for OIF object cell_type = oif.OifCellType( nodes_file=str(tests_common.data_path("sphere393nodes.dat")), triangles_file=str( tests_common.data_path("sphere393triangles.dat")), - system=system, ks=1.0, kb=1.0, kal=1.0, kag=0.1, kv=0.1, + system=self.system, **kwargs, kb=1.0, kal=1.0, kag=0.1, kv=0.1, check_orientation=False, resize=(3.0, 3.0, 3.0)) # creating the OIF object cell0 = oif.OifCell( - cell_type=cell_type, particle_type=0, origin=[5.0, 5.0, 5.0]) - self.assertEqual(system.max_oif_objects, 1) - partcls = system.part.all() + cell_type=cell_type, particle_type=0, origin=half_box_l) + partcls = self.system.part.all() - # fluid diameter_init = cell0.diameter() - print(f"initial diameter = {diameter_init}") + self.assertAlmostEqual(diameter_init, 24., delta=1e-3) # OIF object is being stretched by factor 1.5 - partcls.pos = (partcls.pos - 5) * 1.5 + 5 + partcls.pos = (partcls.pos - half_box_l) * 1.5 + half_box_l diameter_stretched = cell0.diameter() - print(f"stretched diameter = {diameter_stretched}") + self.assertAlmostEqual(diameter_stretched, 36., delta=1e-3) # Apply non-isotropic deformation partcls.pos = partcls.pos * np.array((0.96, 1.05, 1.02)) # Test that restoring forces net to zero and don't produce a torque - system.integrator.run(1) + self.system.integrator.run(1) total_force = np.sum(partcls.f, axis=0) - np.testing.assert_allclose(total_force, [0., 0., 0.], atol=1E-12) + np.testing.assert_allclose(total_force, [0., 0., 0.], atol=1E-10) total_torque = np.zeros(3) - for p in system.part: + for p in self.system.part: total_torque += np.cross(p.pos, p.f) - np.testing.assert_allclose(total_torque, [0., 0., 0.], atol=2E-12) - - # main integration loop - system.thermostat.set_langevin(kT=0, gamma=0.7, seed=42) - # OIF object is let to relax into relaxed shape of the sphere - for _ in range(2): - system.integrator.run(steps=240) - diameter_final = cell0.diameter() - print(f"final diameter = {diameter_final}") - self.assertAlmostEqual( - diameter_final / diameter_init - 1, 0, delta=0.005) + np.testing.assert_allclose(total_torque, [0., 0., 0.], atol=1E-10) + + # Test relaxation to equilibrium diameter (overdamped system) + self.system.thermostat.set_langevin(kT=0, gamma=0.7, seed=42) + # warmup + self.system.integrator.run(steps=50) + # sampling + xdata = [] + ydata = [] + for _ in range(20): + self.system.integrator.run(steps=20) + xdata.append(self.system.time) + ydata.append(cell0.diameter()) + # check exponential decay + (prefactor, lam, _, diameter_final), _ = scipy.optimize.curve_fit( + lambda x, a, b, c, d: a * np.exp(-b * x + c) + d, xdata, ydata, + p0=[3., 0.03, 0., diameter_init], + bounds=([-np.inf, 0., -np.inf, 0.], 4 * [np.inf])) + self.assertGreater(prefactor, 0.) + self.assertAlmostEqual(diameter_final, diameter_init, delta=0.005) + self.assertAlmostEqual(lam, 0.0325, delta=0.0001) + + def test(self): + self.assertEqual(self.system.max_oif_objects, 0) + with self.subTest(msg='linear stretching'): + self.check_relaxation(kslin=1.) + self.assertEqual(self.system.max_oif_objects, 1) + with self.subTest(msg='neo-Hookean stretching'): + self.check_relaxation(ks=1.) + self.assertEqual(self.system.max_oif_objects, 2) if __name__ == "__main__": From 34d9d42a82153504189038d694163e01c7722744 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 15 Dec 2022 21:14:27 +0100 Subject: [PATCH 5/6] core: Remove unreachable code --- src/core/MpiCallbacks.hpp | 6 ++-- src/core/p3m/common.cpp | 51 +++++++++++-------------------- src/utils/include/utils/index.hpp | 8 ++--- 3 files changed, 22 insertions(+), 43 deletions(-) diff --git a/src/core/MpiCallbacks.hpp b/src/core/MpiCallbacks.hpp index 78a27a5b1d0..7722cdee5de 100644 --- a/src/core/MpiCallbacks.hpp +++ b/src/core/MpiCallbacks.hpp @@ -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); diff --git a/src/core/p3m/common.cpp b/src/core/p3m/common.cpp index 7d010576bd8..bf5bdf4f4b4 100644 --- a/src/core/p3m/common.cpp +++ b/src/core/p3m/common.cpp @@ -26,68 +26,53 @@ #include "common.hpp" #include "LocalBox.hpp" -#include "communication.hpp" -#include "errorhandling.hpp" #include #include #include #include -#include +#include double p3m_analytic_cotangent_sum(int n, double mesh_i, int cao) { auto const c = Utils::sqr(std::cos(Utils::pi() * mesh_i * static_cast(n))); - auto res = 0.0; switch (cao) { case 1: { - res = 1.0; - break; + return 1.0; } case 2: { - res = (1.0 + c * 2.0) / 3.0; - break; + return (1.0 + c * 2.0) / 3.0; } case 3: { - res = (2.0 + c * (11.0 + c * 2.0)) / 15.0; - break; + return (2.0 + c * (11.0 + c * 2.0)) / 15.0; } case 4: { - res = (17.0 + c * (180.0 + c * (114.0 + c * 4.0))) / 315.0; - break; + return (17.0 + c * (180.0 + c * (114.0 + c * 4.0))) / 315.0; } case 5: { - res = (62.0 + c * (1072.0 + c * (1452.0 + c * (247.0 + c * 2.0)))) / 2835.0; - break; + return (62.0 + c * (1072.0 + c * (1452.0 + c * (247.0 + c * 2.0)))) / + 2835.0; } case 6: { - res = (1382.0 + - c * (35396.0 + - c * (83021.0 + c * (34096.0 + c * (2026.0 + c * 4.0))))) / - 155925.0; - break; + return (1382.0 + + c * (35396.0 + + c * (83021.0 + c * (34096.0 + c * (2026.0 + c * 4.0))))) / + 155925.0; } case 7: { - res = - (21844.0 + - c * (776661.0 + c * (2801040.0 + - c * (2123860.0 + - c * (349500.0 + c * (8166.0 + c * 4.0)))))) / - 6081075.0; - break; + return (21844.0 + + c * (776661.0 + + c * (2801040.0 + + c * (2123860.0 + + c * (349500.0 + c * (8166.0 + c * 4.0)))))) / + 6081075.0; } default: { - fprintf(stderr, - "%d: INTERNAL_ERROR: The value %d for the interpolation order " - "should not occur!\n", - this_node, cao); - errexit(); + throw std::logic_error("Invalid value cao=" + std::to_string(cao)); } } - - return res; } void P3MLocalMesh::calc_local_ca_mesh(P3MParameters const ¶ms, diff --git a/src/utils/include/utils/index.hpp b/src/utils/include/utils/index.hpp index bd3930725e5..68b9bbd3f8e 100644 --- a/src/utils/include/utils/index.hpp +++ b/src/utils/include/utils/index.hpp @@ -46,14 +46,10 @@ get_linear_index(int a, int b, int c, const Vector3i &adim, assert((b >= 0) && (b < adim[1])); assert((c >= 0) && (c < adim[2])); - switch (memory_order) { - case MemoryOrder::COLUMN_MAJOR: + if (memory_order == MemoryOrder::COLUMN_MAJOR) { return a + adim[0] * (b + adim[1] * c); - case MemoryOrder::ROW_MAJOR: - return adim[1] * adim[2] * a + adim[2] * b + c; - default: - throw std::runtime_error("Unknown memory order"); } + return adim[1] * adim[2] * a + adim[2] * b + c; } inline int From 4ac053765c5540c17b74e564bc4656ca82ec20f9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 16 Dec 2022 01:20:32 +0100 Subject: [PATCH 6/6] Improve code coverage Factor out duplicated code, remove unreachable code, add extra checks. --- .gitlab-ci.yml | 1 + src/core/MpiCallbacks.hpp | 2 +- src/core/unit_tests/MpiCallbacks_test.cpp | 1 + .../Particle_serialization_test.cpp | 4 +- src/core/unit_tests/Verlet_list_test.cpp | 7 +- src/core/unit_tests/p3m_test.cpp | 26 ++++++ src/utils/tests/Factory_test.cpp | 49 +++------- src/utils/tests/gatherv_test.cpp | 89 +++++++------------ testsuite/python/ek_eof_one_species.py | 51 +++-------- testsuite/python/integrator_npt.py | 43 ++++----- 10 files changed, 116 insertions(+), 157 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index d7b0a2a7efa..c0789c2a4c2 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -292,6 +292,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' diff --git a/src/core/MpiCallbacks.hpp b/src/core/MpiCallbacks.hpp index 7722cdee5de..03f3a014b79 100644 --- a/src/core/MpiCallbacks.hpp +++ b/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. */ diff --git a/src/core/unit_tests/MpiCallbacks_test.cpp b/src/core/unit_tests/MpiCallbacks_test.cpp index 5f70428981c..3afbe3f2a63 100644 --- a/src/core/unit_tests/MpiCallbacks_test.cpp +++ b/src/core/unit_tests/MpiCallbacks_test.cpp @@ -33,6 +33,7 @@ #include #include +#include #include static bool called = false; diff --git a/src/core/unit_tests/Particle_serialization_test.cpp b/src/core/unit_tests/Particle_serialization_test.cpp index 43da852db18..0e6ebec1658 100644 --- a/src/core/unit_tests/Particle_serialization_test.cpp +++ b/src/core/unit_tests/Particle_serialization_test.cpp @@ -213,7 +213,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE( { static_assert(!TraitWrapper::template apply::value); typename Checker::buffer_type buffer_ref = {"Particle"}; - typename Checker::buffer_type buffer; + typename Checker::buffer_type buffer = {}; Checker oa{buffer}; Particle p; oa &p; @@ -228,7 +228,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE( "Utils::compact_vector", #endif }; - typename Checker::buffer_type buffer; + typename Checker::buffer_type buffer = {}; Checker oa{buffer}; Particle p; oa | p; diff --git a/src/core/unit_tests/Verlet_list_test.cpp b/src/core/unit_tests/Verlet_list_test.cpp index a83488a7a45..9749a70f42d 100644 --- a/src/core/unit_tests/Verlet_list_test.cpp +++ b/src/core/unit_tests/Verlet_list_test.cpp @@ -280,8 +280,11 @@ int main(int argc, char **argv) { espresso::system = std::make_unique(argc, argv); // the test case only works for 4 MPI ranks boost::mpi::communicator world; - if (world.size() == 4) - return boost::unit_test::unit_test_main(init_unit_test, argc, argv); + int error_code = 0; + if (world.size() == 4) { + error_code = boost::unit_test::unit_test_main(init_unit_test, argc, argv); + } + return error_code; } #else // ifdef LENNARD_JONES int main(int argc, char **argv) {} diff --git a/src/core/unit_tests/p3m_test.cpp b/src/core/unit_tests/p3m_test.cpp index 27c88487123..8279ccc6e2e 100644 --- a/src/core/unit_tests/p3m_test.cpp +++ b/src/core/unit_tests/p3m_test.cpp @@ -27,6 +27,8 @@ #include #include +#include +#include #include BOOST_AUTO_TEST_CASE(calc_meshift_false) { @@ -58,3 +60,27 @@ BOOST_AUTO_TEST_CASE(calc_meshift_true) { } } } + +#if defined(P3M) || defined(DP3M) +BOOST_AUTO_TEST_CASE(analytic_cotangent_sum) { + auto constexpr kernel = p3m_analytic_cotangent_sum; + auto constexpr tol = 100. * std::numeric_limits::epsilon(); + + // check only trivial cases + for (auto const cao : {1, 2, 3, 4, 5, 6, 7}) { + BOOST_CHECK_CLOSE(kernel(0, 0., cao), 1., tol); + } + BOOST_CHECK_CLOSE(kernel(1, 0.5, 1), 1., tol); + BOOST_CHECK_CLOSE(kernel(1, 0.5, 2), 1. / 3., tol); + BOOST_CHECK_CLOSE(kernel(1, 0.5, 3), 2. / 15., tol); + BOOST_CHECK_CLOSE(kernel(1, 0.5, 4), 17. / 315., tol); + BOOST_CHECK_CLOSE(kernel(1, 0.5, 5), 62. / 2835., tol); + BOOST_CHECK_CLOSE(kernel(1, 0.5, 6), 1382. / 155925., tol); + BOOST_CHECK_CLOSE(kernel(1, 0.5, 7), 21844. / 6081075., tol); + + // check assertion + for (auto const invalid_cao : {-1, 0, 8}) { + BOOST_CHECK_THROW(kernel(1, 0., invalid_cao), std::logic_error); + } +} +#endif // defined(P3M) || defined(DP3M) diff --git a/src/utils/tests/Factory_test.cpp b/src/utils/tests/Factory_test.cpp index eab88b60a3b..6b853f6c289 100644 --- a/src/utils/tests/Factory_test.cpp +++ b/src/utils/tests/Factory_test.cpp @@ -17,12 +17,6 @@ * along with this program. If not, see . */ -/* Unit tests for the Utils::Factory class. - * The factory is tested by registering different types of classes - * with it (first test), and then checking if instances of those classes can be - * made via the Factory (second test). - */ - #define BOOST_TEST_MODULE Factory test #define BOOST_TEST_DYN_LINK #include @@ -41,43 +35,28 @@ struct DerivedTestClass : public TestClass { void method() override {} }; -struct OtherDerivedTestClass : public TestClass { - void method() override {} -}; - -/* Check registration of construction functions */ -BOOST_AUTO_TEST_CASE(register_class) { - Utils::Factory factory; - - factory.register_new("other_derived_class"); - - BOOST_CHECK(factory.has_builder("other_derived_class")); -} - -/* Check object construction. */ BOOST_AUTO_TEST_CASE(make) { Utils::Factory factory; - factory.register_new("derived_test_class"); + auto const derived_class_name = std::string{"derived_test_class"}; - /* Make a derived object */ - auto o = factory.make("derived_test_class"); - BOOST_CHECK(o); - - /* Check for correct (derived) type */ - BOOST_CHECK(dynamic_cast(o.get()) != nullptr); -} - -BOOST_AUTO_TEST_CASE(type_name) { - const std::string derived_class_name = "derived_test_class"; - - Utils::Factory factory; + // Register construction function factory.register_new(derived_class_name); - /* Make an object */ + // Check registration of construction function + BOOST_REQUIRE(factory.has_builder(derived_class_name)); + + // Make a derived object auto o = factory.make(derived_class_name); + BOOST_REQUIRE(o); o->method(); + + // Check for correct type name BOOST_CHECK_EQUAL(factory.type_name(*o.get()), derived_class_name); - /* Make an unknown object */ + // Check for correct (derived) type + BOOST_CHECK(dynamic_cast(o.get()) != nullptr); + + // Make an unknown object + BOOST_CHECK(not factory.has_builder("unknown")); BOOST_CHECK_THROW(factory.make("unknown"), std::domain_error); } diff --git a/src/utils/tests/gatherv_test.cpp b/src/utils/tests/gatherv_test.cpp index 96ca48925fb..1cd67ba534f 100644 --- a/src/utils/tests/gatherv_test.cpp +++ b/src/utils/tests/gatherv_test.cpp @@ -28,56 +28,72 @@ #include #include +#include #include -/* - * Check that implementation behaves - * like @c MPI_Gatherv with an mpi datatype. - * To test this, we gather the rank from - * every rank to one rank and then check - * that the value was written to the - * correct position in the output array. - */ -BOOST_AUTO_TEST_CASE(mpi_type) { +struct identity { + template constexpr T &&operator()(T &&t) const noexcept { + return std::forward(t); + } +}; + +template void check(T default_value, F conversion) { boost::mpi::communicator world; auto const rank = world.rank(); auto const size = world.size(); auto const root = world.size() - 1; + auto const in = conversion(rank); /* out-of-place */ { if (rank == root) { - std::vector out(size, -1); + std::vector out(size, default_value); std::vector sizes(size, 1); - Utils::Mpi::gatherv(world, &rank, 1, out.data(), sizes.data(), root); + Utils::Mpi::gatherv(world, &in, 1, out.data(), sizes.data(), root); for (int i = 0; i < size; i++) { - BOOST_CHECK_EQUAL(out.at(i), i); + BOOST_CHECK_EQUAL(out.at(i), conversion(i)); } + } else if (rank % 2 == 0) { + Utils::Mpi::gatherv(world, &in, 1, static_cast(nullptr), nullptr, + root); } else { - Utils::Mpi::gatherv(world, &rank, 1, root); + Utils::Mpi::gatherv(world, &in, 1, root); } } /* in-place */ { if (rank == root) { - std::vector out(size, -1); - out[rank] = rank; + std::vector out(size, default_value); + out[rank] = in; std::vector sizes(size, 1); Utils::Mpi::gatherv(world, out.data(), 1, out.data(), sizes.data(), root); for (int i = 0; i < size; i++) { - BOOST_CHECK_EQUAL(out.at(i), i); + BOOST_CHECK_EQUAL(out.at(i), conversion(i)); } + } else if (rank % 2 == 0) { + Utils::Mpi::gatherv(world, &in, 1, static_cast(nullptr), nullptr, + root); } else { - Utils::Mpi::gatherv(world, &rank, 1, root); + Utils::Mpi::gatherv(world, &in, 1, root); } } } +/* + * Check that implementation behaves + * like @c MPI_Gatherv with an mpi datatype. + * To test this, we gather the rank from + * every rank to one rank and then check + * that the value was written to the + * correct position in the output array. + */ +BOOST_AUTO_TEST_CASE(mpi_type) { check(-1, identity{}); } + /* * Check that implementation behaves * like @c MPI_Gatherv with a non-mpi datatype. @@ -87,44 +103,7 @@ BOOST_AUTO_TEST_CASE(mpi_type) { * correct position in the output array. */ BOOST_AUTO_TEST_CASE(non_mpi_type) { - boost::mpi::communicator world; - auto const rank = world.rank(); - auto const size = world.size(); - auto const root = world.size() - 1; - auto const in = std::to_string(rank); - - /* out-of-place */ - { - if (rank == root) { - std::vector out(size); - std::vector sizes(size, 1); - - Utils::Mpi::gatherv(world, &in, 1, out.data(), sizes.data(), root); - - for (int i = 0; i < size; i++) { - BOOST_CHECK_EQUAL(out.at(i), std::to_string(i)); - } - } else { - Utils::Mpi::gatherv(world, &in, 1, root); - } - } - - /* in-place */ - { - if (rank == root) { - std::vector out(size); - out[rank] = in; - std::vector sizes(size, 1); - - Utils::Mpi::gatherv(world, out.data(), 1, out.data(), sizes.data(), root); - - for (int i = 0; i < size; i++) { - BOOST_CHECK_EQUAL(out.at(i), std::to_string(i)); - } - } else { - Utils::Mpi::gatherv(world, &in, 1, root); - } - } + check(std::string{""}, static_cast(std::to_string)); } int main(int argc, char **argv) { diff --git a/testsuite/python/ek_eof_one_species.py b/testsuite/python/ek_eof_one_species.py index 0d3690facf8..ca7753e0708 100644 --- a/testsuite/python/ek_eof_one_species.py +++ b/testsuite/python/ek_eof_one_species.py @@ -24,7 +24,6 @@ import tempfile import contextlib -import sys import math import numpy as np @@ -159,55 +158,25 @@ def hydrostatic_pressure(ek, tensor_entry, box_x, box_y, box_z, agrid): def bisection(): + args = [params_base[k] + for k in ('width', 'bjerrum_length', 'sigma', 'valency')] # initial parameters for bisection scheme size = math.pi / (2.0 * params_base['width']) pnt0 = 0.0 pntm = pnt0 + size pnt1 = pnt0 + 1.9 * size - # the bisection scheme tol = 1.0e-08 while size > tol: - val0 = solve( - pnt0, - params_base['width'], - params_base['bjerrum_length'], - params_base['sigma'], - params_base['valency']) - val1 = solve( - pnt1, - params_base['width'], - params_base['bjerrum_length'], - params_base['sigma'], - params_base['valency']) - valm = solve( - pntm, - params_base['width'], - params_base['bjerrum_length'], - params_base['sigma'], - params_base['valency']) - - if (val0 < 0.0 and val1 > 0.0): - if valm < 0.0: - pnt0 = pntm - size = size / 2.0 - pntm = pnt0 + size - else: - pnt1 = pntm - size = size / 2.0 - pntm = pnt1 - size - elif (val0 > 0.0 and val1 < 0.0): - if valm < 0.0: - pnt1 = pntm - size = size / 2.0 - pntm = pnt1 - size - else: - pnt0 = pntm - size = size / 2.0 - pntm = pnt0 + size + size /= 2.0 + val0, val1, valm = map(lambda x: solve(x, *args), [pnt0, pnt1, pntm]) + assert val0 < 0.0 and val1 > 0.0, "Bisection method failed" + if valm < 0.0: + pnt0 = pntm + pntm += size else: - sys.exit("Bisection method fails:\n" - "Tuning of regular boundaries may be required.") + pnt1 = pntm + pntm -= size return pntm diff --git a/testsuite/python/integrator_npt.py b/testsuite/python/integrator_npt.py index 3007308a705..8684db9966f 100644 --- a/testsuite/python/integrator_npt.py +++ b/testsuite/python/integrator_npt.py @@ -120,8 +120,11 @@ def test_00_integrator_recovery(self): np.copy(system.part.all().pos), positions_start + np.array([[-1.2e-3, 0, 0], [1.2e-3, 0, 0]])) - def run_with_p3m(self, p3m, **npt_kwargs): + def run_with_p3m(self, p3m, method): system = self.system + npt_kwargs = {"ext_pressure": 0.001, "piston": 0.001} + npt_kwargs_rectangular = { + "cubic_box": False, "direction": (False, True, True), **npt_kwargs} np.random.seed(42) # set up particles system.box_l = [6] * 3 @@ -140,10 +143,23 @@ def run_with_p3m(self, p3m, **npt_kwargs): # combine NpT with a P3M algorithm system.actors.add(p3m) system.integrator.run(10) - system.integrator.set_isotropic_npt(ext_pressure=0.001, piston=0.001, - **npt_kwargs) + system.integrator.set_isotropic_npt(**npt_kwargs) system.thermostat.set_npt(kT=1.0, gamma0=2, gammav=0.04, seed=42) system.integrator.run(10) + # check runtime warnings + self.system.thermostat.turn_off() + self.system.integrator.set_vv() + err_msg = f"If {method} is being used you must use the cubic box NpT" + with self.assertRaisesRegex(RuntimeError, err_msg): + system.integrator.set_isotropic_npt(**npt_kwargs_rectangular) + self.assertIsInstance( + system.integrator.integrator, + espressomd.integrate.VelocityVerlet) + system.actors.remove(p3m) + system.integrator.set_isotropic_npt(**npt_kwargs_rectangular) + system.actors.add(p3m) + with self.assertRaisesRegex(Exception, err_msg): + system.integrator.run(0, recalc_forces=True) @utx.skipIfMissingFeatures(["DP3M"]) def test_npt_dp3m_cpu(self): @@ -151,12 +167,7 @@ def test_npt_dp3m_cpu(self): dp3m = espressomd.magnetostatics.DipolarP3M( prefactor=1.0, accuracy=1e-2, mesh=3 * [36], cao=7, r_cut=1.0, alpha=2.995, tune=False) - with self.assertRaisesRegex(RuntimeError, 'If magnetostatics is being used you must use the cubic box NpT'): - self.run_with_p3m( - dp3m, cubic_box=False, direction=(False, True, True)) - self.tearDown() - # should not raise an exception - self.run_with_p3m(dp3m) + self.run_with_p3m(dp3m, "magnetostatics") @utx.skipIfMissingFeatures(["P3M"]) def test_npt_p3m_cpu(self): @@ -164,12 +175,7 @@ def test_npt_p3m_cpu(self): p3m = espressomd.electrostatics.P3M( prefactor=1.0, accuracy=1e-2, mesh=3 * [8], cao=3, r_cut=0.36, alpha=5.35, tune=False) - with self.assertRaisesRegex(RuntimeError, 'If electrostatics is being used you must use the cubic box NpT'): - self.run_with_p3m( - p3m, cubic_box=False, direction=(False, True, True)) - self.tearDown() - # should not raise an exception - self.run_with_p3m(p3m) + self.run_with_p3m(p3m, "electrostatics") @utx.skipIfMissingGPU() @utx.skipIfMissingFeatures(["P3M"]) @@ -178,12 +184,7 @@ def test_npt_p3m_gpu(self): p3m = espressomd.electrostatics.P3MGPU( prefactor=1.0, accuracy=1e-2, mesh=3 * [8], cao=3, r_cut=0.36, alpha=5.35, tune=False) - with self.assertRaisesRegex(RuntimeError, 'If electrostatics is being used you must use the cubic box NpT'): - self.run_with_p3m( - p3m, cubic_box=False, direction=(False, True, True)) - self.tearDown() - # should not raise an exception - self.run_with_p3m(p3m) + self.run_with_p3m(p3m, "electrostatics") if __name__ == "__main__":