From 9771adbfef7e8f017b5af0c5b23ce299632b2b8c Mon Sep 17 00:00:00 2001 From: "kodiakhq[bot]" <49736102+kodiakhq[bot]@users.noreply.github.com> Date: Mon, 17 Apr 2023 15:25:26 +0000 Subject: [PATCH] Add support for ARM architectures (#4708) Description of changes: - adjust test tolerances for architectures that don't have double extended precision, e.g. ARM - add sanity checks in chain analysis functions and handle the corner case `number_of_chains=1` - re-enable support for fast-math mode --- maintainer/CI/build_cmake.sh | 2 +- src/core/electrostatics/icc.cpp | 8 +++- .../tests/ConstantpHEnsemble_test.cpp | 4 +- .../tests/ReactionAlgorithm_test.cpp | 2 +- .../tests/ReactionEnsemble_test.cpp | 4 +- .../tests/SingleReaction_test.cpp | 4 +- .../tests/reaction_methods_utils_test.cpp | 4 +- src/core/statistics_chain.cpp | 31 ++++++------- src/core/statistics_chain.hpp | 15 +++---- .../EspressoSystemStandAlone_test.cpp | 10 ++--- src/core/unit_tests/Verlet_list_test.cpp | 2 +- .../field_coupling_couplings_test.cpp | 3 +- .../unit_tests/field_coupling_fields_test.cpp | 8 +++- .../field_coupling_force_field_test.cpp | 3 +- src/core/unit_tests/p3m_test.cpp | 2 +- src/core/unit_tests/periodic_fold_test.cpp | 4 +- src/core/unit_tests/specfunc_test.cpp | 2 +- src/core/unit_tests/thermostats_test.cpp | 4 +- src/particle_observables/tests/algorithms.cpp | 2 +- src/python/espressomd/analyze.pyx | 14 +++--- src/shapes/unit_tests/Ellipsoid_test.cpp | 2 +- src/shapes/unit_tests/Sphere_test.cpp | 2 +- src/utils/tests/Vector_test.cpp | 5 ++- src/utils/tests/interpolation_test.cpp | 3 +- src/utils/tests/matrix_vector_product.cpp | 5 ++- src/utils/tests/vec_rotate_test.cpp | 2 +- testsuite/python/analyze_chains.py | 43 ++++++++++++++----- testsuite/python/integrator_npt_stats.py | 2 +- 28 files changed, 113 insertions(+), 79 deletions(-) diff --git a/maintainer/CI/build_cmake.sh b/maintainer/CI/build_cmake.sh index 55a34d4d27d..456aeddabe5 100755 --- a/maintainer/CI/build_cmake.sh +++ b/maintainer/CI/build_cmake.sh @@ -123,7 +123,7 @@ if [ "${with_coverage}" = true ]; then fi if [ "${with_fast_math}" = true ]; then - cmake_param_protected="-DCMAKE_CXX_FLAGS=-ffast-math -fno-finite-math-only" + cmake_param_protected="-DCMAKE_CXX_FLAGS=-ffast-math" fi cmake_params="-DCMAKE_BUILD_TYPE=${build_type} -DCMAKE_CXX_STANDARD=${with_cxx_standard} -DWARNINGS_ARE_ERRORS=ON ${cmake_params}" diff --git a/src/core/electrostatics/icc.cpp b/src/core/electrostatics/icc.cpp index 4f555c4841e..e3883dba609 100644 --- a/src/core/electrostatics/icc.cpp +++ b/src/core/electrostatics/icc.cpp @@ -129,6 +129,12 @@ void ICCStar::iteration(CellStructure &cell_structure, for (auto &p : particles) { auto const pid = p.id(); if (pid >= icc_cfg.first_id and pid < icc_cfg.n_icc + icc_cfg.first_id) { + if (p.q() == 0.) { + runtimeErrorMsg() + << "ICC found zero electric charge on a particle. This must " + "never happen"; + break; + } auto const id = p.id() - icc_cfg.first_id; /* the dielectric-related prefactor: */ auto const eps_in = icc_cfg.epsilons[id]; @@ -176,7 +182,7 @@ void ICCStar::iteration(CellStructure &cell_structure, << "Particle with id " << p.id() << " has a charge (q=" << p.q() << ") that is too large for the ICC algorithm"; - max_rel_diff = std::numeric_limits::infinity(); + max_rel_diff = std::numeric_limits::max(); break; } } diff --git a/src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp b/src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp index 126cf851f7e..ddb42ab33a0 100644 --- a/src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp +++ b/src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp @@ -47,7 +47,7 @@ BOOST_AUTO_TEST_CASE(ConstantpHEnsemble_test) { using ConstantpHEnsemble::calculate_acceptance_probability; using ConstantpHEnsemble::ConstantpHEnsemble; }; - constexpr double tol = 100 * std::numeric_limits::epsilon(); + auto constexpr tol = 8. * 100. * std::numeric_limits::epsilon(); ConstantpHEnsembleTest r_algo(42, 20., 0., 1., {}); @@ -78,7 +78,7 @@ BOOST_AUTO_TEST_CASE(ConstantpHEnsemble_test) { std::log10(reaction.gamma))); auto const acceptance = r_algo.calculate_acceptance_probability( reaction, energy, 0., p_numbers); - BOOST_CHECK_CLOSE(acceptance, acceptance_ref, 5 * tol); + BOOST_CHECK_CLOSE(acceptance, acceptance_ref, 5. * tol); } } } diff --git a/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp b/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp index 265746cd7d6..952bcbe9150 100644 --- a/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp +++ b/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp @@ -59,7 +59,7 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { using ReactionAlgorithm::get_random_position_in_box; using ReactionAlgorithm::ReactionAlgorithm; }; - constexpr double tol = 100 * std::numeric_limits::epsilon(); + auto constexpr tol = 8. * 100. * std::numeric_limits::epsilon(); // check acceptance rate ReactionAlgorithmTest r_algo(42, 1., 0., {}); diff --git a/src/core/reaction_methods/tests/ReactionEnsemble_test.cpp b/src/core/reaction_methods/tests/ReactionEnsemble_test.cpp index 1c6934a87ba..e8fbc9decc5 100644 --- a/src/core/reaction_methods/tests/ReactionEnsemble_test.cpp +++ b/src/core/reaction_methods/tests/ReactionEnsemble_test.cpp @@ -57,7 +57,7 @@ BOOST_AUTO_TEST_CASE(ReactionEnsemble_test) { using ReactionEnsemble::generic_oneway_reaction; using ReactionEnsemble::ReactionEnsemble; }; - auto constexpr tol = 100 * std::numeric_limits::epsilon(); + auto constexpr tol = 8. * 100. * std::numeric_limits::epsilon(); // check basic interface { @@ -89,7 +89,7 @@ BOOST_AUTO_TEST_CASE(ReactionEnsemble_test) { std::exp(energy / r_algo.kT); auto const acceptance = r_algo.calculate_acceptance_probability( reaction, energy, 0., p_numbers); - BOOST_CHECK_CLOSE(acceptance, acceptance_ref, 5 * tol); + BOOST_CHECK_CLOSE(acceptance, acceptance_ref, 5. * tol); } } } diff --git a/src/core/reaction_methods/tests/SingleReaction_test.cpp b/src/core/reaction_methods/tests/SingleReaction_test.cpp index 6d9c6a43e5b..0e074c8f23f 100644 --- a/src/core/reaction_methods/tests/SingleReaction_test.cpp +++ b/src/core/reaction_methods/tests/SingleReaction_test.cpp @@ -34,7 +34,7 @@ // and the configurational move probability for a given system state. BOOST_AUTO_TEST_CASE(SingleReaction_test) { using namespace ReactionMethods; - constexpr double tol = 100 * std::numeric_limits::epsilon(); + auto constexpr tol = 8. * 100. * std::numeric_limits::epsilon(); // create a reaction A -> 3 B + 4 C int const type_A = 0; @@ -66,7 +66,7 @@ BOOST_AUTO_TEST_CASE(SingleReaction_test) { std::map{{type_A, i}, {type_B, j}, {type_C, k}}; auto const val = calculate_factorial_expression(reaction, p_numbers); auto const ref = g(i, -1) * g(j, 3) * g(k, 4); - BOOST_CHECK_CLOSE(val, ref, 5 * tol); + BOOST_CHECK_CLOSE(val, ref, 5. * tol); } } } diff --git a/src/core/reaction_methods/tests/reaction_methods_utils_test.cpp b/src/core/reaction_methods/tests/reaction_methods_utils_test.cpp index 579fd206d34..bd9d8b5db34 100644 --- a/src/core/reaction_methods/tests/reaction_methods_utils_test.cpp +++ b/src/core/reaction_methods/tests/reaction_methods_utils_test.cpp @@ -32,7 +32,7 @@ BOOST_AUTO_TEST_CASE(factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i_test) { using namespace ReactionMethods; - constexpr double tol = 100 * std::numeric_limits::epsilon(); + auto constexpr tol = 8. * 100. * std::numeric_limits::epsilon(); auto const reaction_ensemble_combinations = [](int N, int nu) { return (N + nu < 0) ? 0. : std::tgamma(N + 1) / std::tgamma(N + nu + 1); @@ -42,7 +42,7 @@ BOOST_AUTO_TEST_CASE(factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i_test) { for (int nu = -4; nu <= 4; ++nu) { auto const val = factorial_Ni0_divided_by_factorial_Ni0_plus_nu_i(N0, nu); auto const ref = reaction_ensemble_combinations(N0, nu); - BOOST_CHECK_CLOSE(val, ref, 10 * tol); + BOOST_CHECK_CLOSE(val, ref, 10. * tol); } } } diff --git a/src/core/statistics_chain.cpp b/src/core/statistics_chain.cpp index f5122bd005d..575a5784f0d 100644 --- a/src/core/statistics_chain.cpp +++ b/src/core/statistics_chain.cpp @@ -34,12 +34,11 @@ #include #include -std::array calc_re(int chain_start, int chain_n_chains, - int chain_length) { +std::array calc_re(int chain_start, int n_chains, int chain_length) { double dist = 0.0, dist2 = 0.0, dist4 = 0.0; std::array re; - for (int i = 0; i < chain_n_chains; i++) { + for (int i = 0; i < n_chains; i++) { auto const &p1 = get_particle_data(chain_start + i * chain_length + chain_length - 1); auto const &p2 = get_particle_data(chain_start + i * chain_length); @@ -52,20 +51,19 @@ std::array calc_re(int chain_start, int chain_n_chains, dist2 += norm2; dist4 += norm2 * norm2; } - auto const tmp = static_cast(chain_n_chains); + auto const tmp = static_cast(n_chains); re[0] = dist / tmp; re[2] = dist2 / tmp; - re[1] = sqrt(re[2] - re[0] * re[0]); - re[3] = sqrt(dist4 / tmp - re[2] * re[2]); + re[1] = (n_chains == 1) ? 0. : std::sqrt(re[2] - Utils::sqr(re[0])); + re[3] = (n_chains == 1) ? 0. : std::sqrt(dist4 / tmp - Utils::sqr(re[2])); return re; } -std::array calc_rg(int chain_start, int chain_n_chains, - int chain_length) { +std::array calc_rg(int chain_start, int n_chains, int chain_length) { double r_G = 0.0, r_G2 = 0.0, r_G4 = 0.0; std::array rg; - for (int i = 0; i < chain_n_chains; i++) { + for (int i = 0; i < n_chains; i++) { double M = 0.0; Utils::Vector3d r_CM{}; for (int j = 0; j < chain_length; j++) { @@ -93,22 +91,21 @@ std::array calc_rg(int chain_start, int chain_n_chains, r_G2 += tmp; r_G4 += tmp * tmp; } - auto const tmp = static_cast(chain_n_chains); + auto const tmp = static_cast(n_chains); rg[0] = r_G / tmp; rg[2] = r_G2 / tmp; - rg[1] = sqrt(rg[2] - Utils::sqr(rg[0])); - rg[3] = sqrt(r_G4 / tmp - Utils::sqr(rg[2])); + rg[1] = (n_chains == 1) ? 0. : std::sqrt(rg[2] - Utils::sqr(rg[0])); + rg[3] = (n_chains == 1) ? 0. : std::sqrt(r_G4 / tmp - Utils::sqr(rg[2])); return rg; } -std::array calc_rh(int chain_start, int chain_n_chains, - int chain_length) { +std::array calc_rh(int chain_start, int n_chains, int chain_length) { double r_H = 0.0, r_H2 = 0.0; std::array rh; auto const chain_l = static_cast(chain_length); auto const prefac = 0.5 * chain_l * (chain_l - 1.); - for (int p = 0; p < chain_n_chains; p++) { + for (int p = 0; p < n_chains; p++) { double ri = 0.0; for (int i = chain_start + chain_length * p; i < chain_start + chain_length * (p + 1); i++) { @@ -125,8 +122,8 @@ std::array calc_rh(int chain_start, int chain_n_chains, r_H += tmp; r_H2 += tmp * tmp; } - auto const tmp = static_cast(chain_n_chains); + auto const tmp = static_cast(n_chains); rh[0] = r_H / tmp; - rh[1] = sqrt(r_H2 / tmp - Utils::sqr(rh[0])); + rh[1] = (n_chains == 1) ? 0. : std::sqrt(r_H2 / tmp - Utils::sqr(rh[0])); return rh; } diff --git a/src/core/statistics_chain.hpp b/src/core/statistics_chain.hpp index 63647ae711d..8d420b6b6bb 100644 --- a/src/core/statistics_chain.hpp +++ b/src/core/statistics_chain.hpp @@ -38,11 +38,10 @@ * of monodisperse polymers with continuous ids. * * @param chain_start The id of the first monomer of the first chain. - * @param chain_n_chains Number of chains contained in the range. + * @param n_chains Number of chains contained in the range. * @param chain_length The length of every chain. */ -std::array calc_re(int chain_start, int chain_n_chains, - int chain_length); +std::array calc_re(int chain_start, int n_chains, int chain_length); /** * @brief Calculate the radius of gyration. @@ -51,11 +50,10 @@ std::array calc_re(int chain_start, int chain_n_chains, * of monodisperse polymers with continuous ids. * * @param chain_start The id of the first monomer of the first chain. - * @param chain_n_chains Number of chains contained in the range. + * @param n_chains Number of chains contained in the range. * @param chain_length The length of every chain. */ -std::array calc_rg(int chain_start, int chain_n_chains, - int chain_length); +std::array calc_rg(int chain_start, int n_chains, int chain_length); /** * @brief Calculate the hydrodynamic radius (ref. Kirkwood-Zimm theory). @@ -64,11 +62,10 @@ std::array calc_rg(int chain_start, int chain_n_chains, * of monodisperse polymers with continuous ids. * * @param chain_start The id of the first monomer of the first chain. - * @param chain_n_chains Number of chains contained in the range. + * @param n_chains Number of chains contained in the range. * @param chain_length The length of every chain. */ -std::array calc_rh(int chain_start, int chain_n_chains, - int chain_length); +std::array calc_rh(int chain_start, int n_chains, int chain_length); /**@}*/ #endif diff --git a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp index f2cdabff335..1dab14defe7 100644 --- a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp +++ b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp @@ -118,7 +118,7 @@ static void mpi_set_tuned_p3m(double prefactor) { BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory, *utf::precondition(if_head_node())) { - constexpr auto tol = 100. * std::numeric_limits::epsilon(); + auto constexpr tol = 8. * 100. * std::numeric_limits::epsilon(); auto const box_l = 8.; auto const box_center = box_l / 2.; @@ -226,8 +226,8 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory, for (int i = 0; i < n_pairs; ++i) { auto const ref_inter = (i == lj_pair_ab) ? lj_energy : 0.; auto const ref_intra = (i == lj_pair_bb) ? lj_energy : 0.; - BOOST_CHECK_CLOSE(obs_energy->non_bonded_inter[i], ref_inter, 500. * tol); - BOOST_CHECK_CLOSE(obs_energy->non_bonded_intra[i], ref_intra, 500. * tol); + BOOST_CHECK_CLOSE(obs_energy->non_bonded_inter[i], ref_inter, 1e-10); + BOOST_CHECK_CLOSE(obs_energy->non_bonded_intra[i], ref_intra, 1e-10); } } #endif // LENNARD_JONES @@ -256,8 +256,8 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory, -0.5 * fene_bond.k * Utils::sqr(fene_bond.drmax) * std::log(1.0 - Utils::sqr((dist - fene_bond.r0) / fene_bond.drmax)); BOOST_CHECK_CLOSE(obs_energy->bonded[none_bond_id], none_energy, 0.0); - BOOST_CHECK_CLOSE(obs_energy->bonded[harm_bond_id], harm_energy, 40. * tol); - BOOST_CHECK_CLOSE(obs_energy->bonded[fene_bond_id], fene_energy, 40. * tol); + BOOST_CHECK_CLOSE(obs_energy->bonded[harm_bond_id], harm_energy, 1e-10); + BOOST_CHECK_CLOSE(obs_energy->bonded[fene_bond_id], fene_energy, 1e-10); } // check electrostatics diff --git a/src/core/unit_tests/Verlet_list_test.cpp b/src/core/unit_tests/Verlet_list_test.cpp index a614d800bc6..43368396528 100644 --- a/src/core/unit_tests/Verlet_list_test.cpp +++ b/src/core/unit_tests/Verlet_list_test.cpp @@ -144,7 +144,7 @@ BOOST_TEST_DECORATOR(*utf::precondition(if_head_node())) BOOST_DATA_TEST_CASE_F(ParticleFactory, verlet_list_update, bdata::make(node_grids) * bdata::make(propagators), node_grid, integration_helper) { - constexpr auto tol = 100. * std::numeric_limits::epsilon(); + auto constexpr tol = 8. * 100. * std::numeric_limits::epsilon(); boost::mpi::communicator world; auto const box_l = 8.; diff --git a/src/core/unit_tests/field_coupling_couplings_test.cpp b/src/core/unit_tests/field_coupling_couplings_test.cpp index 26f52e44ef3..94853f80add 100644 --- a/src/core/unit_tests/field_coupling_couplings_test.cpp +++ b/src/core/unit_tests/field_coupling_couplings_test.cpp @@ -16,7 +16,8 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ -#define BOOST_TEST_MODULE AutoParameter test + +#define BOOST_TEST_MODULE Field coupling test for fields #define BOOST_TEST_DYN_LINK #include diff --git a/src/core/unit_tests/field_coupling_fields_test.cpp b/src/core/unit_tests/field_coupling_fields_test.cpp index 91f2e5f3c7b..86498047ec1 100644 --- a/src/core/unit_tests/field_coupling_fields_test.cpp +++ b/src/core/unit_tests/field_coupling_fields_test.cpp @@ -17,7 +17,7 @@ * along with this program. If not, see . */ -#define BOOST_TEST_MODULE AutoParameter test +#define BOOST_TEST_MODULE Field coupling test for fields #define BOOST_TEST_DYN_LINK #include @@ -37,7 +37,7 @@ #include #include -constexpr auto eps = 10 * std::numeric_limits::epsilon(); +auto constexpr eps = 10. * std::numeric_limits::epsilon(); using namespace FieldCoupling::Fields; @@ -286,6 +286,7 @@ BOOST_AUTO_TEST_CASE(interpolated_scalar_field) { BOOST_CHECK_SMALL(std::abs(interpolated_value - field_value), eps); } +#ifndef __FAST_MATH__ /* jacobian value */ { using Utils::Interpolation::bspline_3d_gradient_accumulate; @@ -313,6 +314,7 @@ BOOST_AUTO_TEST_CASE(interpolated_scalar_field) { BOOST_CHECK_SMALL((interpolated_value - field_value).norm(), eps); } +#endif // __FAST_MATH__ } BOOST_AUTO_TEST_CASE(interpolated_vector_field) { @@ -362,6 +364,7 @@ BOOST_AUTO_TEST_CASE(interpolated_vector_field) { BOOST_CHECK_SMALL((interpolated_value - field_value).norm(), eps); } +#ifndef __FAST_MATH__ /* jacobian value */ { using Utils::Interpolation::bspline_3d_gradient_accumulate; @@ -400,4 +403,5 @@ BOOST_AUTO_TEST_CASE(interpolated_vector_field) { BOOST_CHECK_SMALL( (interpolated_value.row<1>() - field_value.row<1>()).norm(), eps); } +#endif // __FAST_MATH__ } diff --git a/src/core/unit_tests/field_coupling_force_field_test.cpp b/src/core/unit_tests/field_coupling_force_field_test.cpp index 9dd11e00cc8..3ec4c188bb9 100644 --- a/src/core/unit_tests/field_coupling_force_field_test.cpp +++ b/src/core/unit_tests/field_coupling_force_field_test.cpp @@ -16,7 +16,8 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ -#define BOOST_TEST_MODULE test + +#define BOOST_TEST_MODULE Field coupling test for force fields #define BOOST_TEST_DYN_LINK #include diff --git a/src/core/unit_tests/p3m_test.cpp b/src/core/unit_tests/p3m_test.cpp index 8279ccc6e2e..0bc343de094 100644 --- a/src/core/unit_tests/p3m_test.cpp +++ b/src/core/unit_tests/p3m_test.cpp @@ -64,7 +64,7 @@ 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(); + auto constexpr tol = 8. * 100. * std::numeric_limits::epsilon(); // check only trivial cases for (auto const cao : {1, 2, 3, 4, 5, 6, 7}) { diff --git a/src/core/unit_tests/periodic_fold_test.cpp b/src/core/unit_tests/periodic_fold_test.cpp index 7346b472c0e..49d69777c94 100644 --- a/src/core/unit_tests/periodic_fold_test.cpp +++ b/src/core/unit_tests/periodic_fold_test.cpp @@ -106,7 +106,7 @@ BOOST_AUTO_TEST_CASE(with_image_count) { BOOST_CHECK(res.first >= 0.); BOOST_CHECK(res.first <= box); BOOST_CHECK(std::abs(res.first - x + res.second * box) <= - std::numeric_limits::epsilon()); + 4. * std::numeric_limits::epsilon()); } /* Corner right */ @@ -117,7 +117,7 @@ BOOST_AUTO_TEST_CASE(with_image_count) { BOOST_CHECK(res.first >= 0.); BOOST_CHECK(res.first < box); BOOST_CHECK(std::abs(res.first - x + res.second * box) <= - std::numeric_limits::epsilon()); + 4. * std::numeric_limits::epsilon()); } } diff --git a/src/core/unit_tests/specfunc_test.cpp b/src/core/unit_tests/specfunc_test.cpp index 96700c71bbd..d41a470f085 100644 --- a/src/core/unit_tests/specfunc_test.cpp +++ b/src/core/unit_tests/specfunc_test.cpp @@ -26,7 +26,7 @@ #include #include -constexpr auto eps = 100. * std::numeric_limits::epsilon(); +auto constexpr eps = 8. * 100. * std::numeric_limits::epsilon(); BOOST_AUTO_TEST_CASE(hurwitz_zeta_function) { constexpr auto max_bits = 54.0; diff --git a/src/core/unit_tests/thermostats_test.cpp b/src/core/unit_tests/thermostats_test.cpp index c0cca417065..dc069746b87 100644 --- a/src/core/unit_tests/thermostats_test.cpp +++ b/src/core/unit_tests/thermostats_test.cpp @@ -41,8 +41,8 @@ #include // multiply by 100 because BOOST_CHECK_CLOSE takes a percentage tolerance, -// and by 6 to account for error accumulation in thermostat functions -constexpr auto tol = 6 * 100 * std::numeric_limits::epsilon(); +// and by 8 to account for error accumulation in thermostat functions +auto constexpr tol = 8. * 100. * std::numeric_limits::epsilon(); Particle particle_factory() { Particle p{}; diff --git a/src/particle_observables/tests/algorithms.cpp b/src/particle_observables/tests/algorithms.cpp index 4a1cef824cf..7fabb7c3ef4 100644 --- a/src/particle_observables/tests/algorithms.cpp +++ b/src/particle_observables/tests/algorithms.cpp @@ -79,7 +79,7 @@ BOOST_AUTO_TEST_CASE(algorithms_integer) { } BOOST_AUTO_TEST_CASE(algorithms_double) { - auto constexpr tol = 100 * std::numeric_limits::epsilon(); + auto constexpr tol = 8. * 100. * std::numeric_limits::epsilon(); std::vector const values{1., 2., 3., 4.}; { auto const res = WeightedAverage()(values); diff --git a/src/python/espressomd/analyze.pyx b/src/python/espressomd/analyze.pyx index 4c9d67628f7..de629412f3b 100644 --- a/src/python/espressomd/analyze.pyx +++ b/src/python/espressomd/analyze.pyx @@ -505,15 +505,19 @@ class Analysis: number_of_chains, 1, int, "number_of_chains=int is a required argument") utils.check_type_or_throw_except( chain_length, 1, int, "chain_length=int is a required argument") + if number_of_chains <= 0: + raise ValueError("Chain analysis needs at least 1 chain") + if chain_length <= 0: + raise ValueError("Chain analysis needs at least 1 bead per chain") id_min = chain_start id_max = chain_start + chain_length * number_of_chains for i in range(id_min, id_max): if not self._system.part.exists(i): - raise ValueError(f'particle with id {i} does not exist\n' - f'cannot perform analysis on the range ' - f'chain_start={chain_start}, number_of_chains=' - f'{number_of_chains}, chain_length={chain_length}\n' - f'please provide a contiguous range of particle ids') + raise RuntimeError(f'particle with id {i} does not exist\n' + f'cannot perform analysis on the range ' + f'chain_start={chain_start}, number_of_chains=' + f'{number_of_chains}, chain_length={chain_length}\n' + f'please provide a contiguous range of particle ids') # # Structure factor diff --git a/src/shapes/unit_tests/Ellipsoid_test.cpp b/src/shapes/unit_tests/Ellipsoid_test.cpp index 2f4331a6147..d9d41874bf7 100644 --- a/src/shapes/unit_tests/Ellipsoid_test.cpp +++ b/src/shapes/unit_tests/Ellipsoid_test.cpp @@ -34,7 +34,7 @@ BOOST_AUTO_TEST_CASE(dist_function) { // multiply by 100 because BOOST_REQUIRE_CLOSE takes a percentage tolerance - auto constexpr tol = std::numeric_limits::epsilon() * 100; + auto constexpr tol = 8. * 100. * std::numeric_limits::epsilon(); double const semiaxes[3] = {3.1, 2.2, 1.3}; Shapes::Ellipsoid e; diff --git a/src/shapes/unit_tests/Sphere_test.cpp b/src/shapes/unit_tests/Sphere_test.cpp index e6eb19aa2f7..61fb916ce7f 100644 --- a/src/shapes/unit_tests/Sphere_test.cpp +++ b/src/shapes/unit_tests/Sphere_test.cpp @@ -34,7 +34,7 @@ void check_distance_function(Shapes::Sphere &s) { Utils::Vector3d vec; double dist; // multiply by 100 because BOOST_REQUIRE_CLOSE takes a percentage tolerance - auto const tol = std::numeric_limits::epsilon() * 100; + auto constexpr tol = 8. * 100. * std::numeric_limits::epsilon(); s.rad() = 1.0; pos = {0., 0., 0.}; diff --git a/src/utils/tests/Vector_test.cpp b/src/utils/tests/Vector_test.cpp index eacf25a10ad..2dc8ec1a8b6 100644 --- a/src/utils/tests/Vector_test.cpp +++ b/src/utils/tests/Vector_test.cpp @@ -119,10 +119,11 @@ BOOST_AUTO_TEST_CASE(test_norm2) { } BOOST_AUTO_TEST_CASE(normalize) { - Utils::Vector3d v{1, 2, 3}; + auto constexpr tol = 8. * 100. * std::numeric_limits::epsilon(); + Utils::Vector3d v{1., 2., 3.}; v.normalize(); - BOOST_CHECK((v.norm2() - 1.0) <= std::numeric_limits::epsilon()); + BOOST_CHECK_CLOSE(v.norm2(), 1.0, tol); } BOOST_AUTO_TEST_CASE(comparison_operators) { diff --git a/src/utils/tests/interpolation_test.cpp b/src/utils/tests/interpolation_test.cpp index 701a7456c89..b452ab32e9d 100644 --- a/src/utils/tests/interpolation_test.cpp +++ b/src/utils/tests/interpolation_test.cpp @@ -135,6 +135,7 @@ BOOST_AUTO_TEST_CASE(sum_of_weights_odd) { } BOOST_AUTO_TEST_CASE(nearest_point) { + auto constexpr tol = 8. * 100. * std::numeric_limits::epsilon(); std::array nmp; double weight; auto save_ind = [&nmp, &weight](const std::array &ind, double w) { @@ -145,7 +146,7 @@ BOOST_AUTO_TEST_CASE(nearest_point) { bspline_3d<1>({.1, .2, .3}, save_ind, {0.5, 0.5, 0.5}, {}); BOOST_CHECK((std::array{{0, 0, 1}} == nmp)); - BOOST_CHECK_CLOSE(weight, 1., 100. * std::numeric_limits::epsilon()); + BOOST_CHECK_CLOSE(weight, 1., tol); } BOOST_AUTO_TEST_CASE(interpolation_points_3) { diff --git a/src/utils/tests/matrix_vector_product.cpp b/src/utils/tests/matrix_vector_product.cpp index a13a61371dd..13e6e8c6526 100644 --- a/src/utils/tests/matrix_vector_product.cpp +++ b/src/utils/tests/matrix_vector_product.cpp @@ -32,10 +32,11 @@ extern constexpr std::array, 3> matrix{ {{{1, 2, 9}}, {{8, 41, 6}}, {{31, 15, 99}}}}; BOOST_AUTO_TEST_CASE(inner_product) { + auto constexpr tol = 8. * 100. * std::numeric_limits::epsilon(); const std::array vector{{0.5, 1.25, 3.1}}; auto const result = Utils::matrix_vector_product(vector); for (int i = 0; i < 3; ++i) { - BOOST_CHECK_CLOSE(result[i], boost::inner_product(matrix[i], vector, 0.0), - 100. * std::numeric_limits::epsilon()); + auto const ref = boost::inner_product(matrix[i], vector, 0.0); + BOOST_CHECK_CLOSE(result[i], ref, tol); } } diff --git a/src/utils/tests/vec_rotate_test.cpp b/src/utils/tests/vec_rotate_test.cpp index c1ad10eabb7..88b2af17983 100644 --- a/src/utils/tests/vec_rotate_test.cpp +++ b/src/utils/tests/vec_rotate_test.cpp @@ -44,7 +44,7 @@ BOOST_AUTO_TEST_CASE(rotation) { auto const is = Utils::vec_rotate(k, t, v); auto const rel_diff = (expected - is).norm() / expected.norm(); - BOOST_CHECK(rel_diff < std::numeric_limits::epsilon()); + BOOST_CHECK(rel_diff < 8. * std::numeric_limits::epsilon()); } BOOST_AUTO_TEST_CASE(angle_between) { diff --git a/testsuite/python/analyze_chains.py b/testsuite/python/analyze_chains.py index 46c6c351995..adf75eb0d6c 100644 --- a/testsuite/python/analyze_chains.py +++ b/testsuite/python/analyze_chains.py @@ -72,8 +72,12 @@ def calc_re(self, num_poly, num_mono): dist = self.system.part.by_ids( head_ids).pos - self.system.part.by_ids(tail_ids).pos dist = np.sum(dist**2, axis=-1) - return np.mean(np.sqrt(dist)), np.std( - np.sqrt(dist)), np.mean(dist), np.std(dist) + return ( + np.mean(np.sqrt(dist)), + 0. if num_mono == 1 else np.std(np.sqrt(dist), ddof=0), + np.mean(dist), + 0. if num_mono == 1 else np.std(dist, ddof=0), + ) # python version of the espresso core function, # does not check mirror distances @@ -87,8 +91,12 @@ def calc_rg(self, num_poly, num_mono): rg2.append(np.var(r, axis=0)) rg2 = np.array(rg2) rg2 = np.sum(rg2, axis=1) - return np.mean(np.sqrt(rg2)), np.std( - np.sqrt(rg2)), np.mean(rg2), np.std(rg2) + return ( + np.mean(np.sqrt(rg2)), + 0. if num_mono == 1 else np.std(np.sqrt(rg2), ddof=0), + np.mean(rg2), + 0. if num_mono == 1 else np.std(rg2, ddof=0), + ) # python version of the espresso core function, # does not check mirror distances @@ -105,7 +113,7 @@ def calc_rh(self, num_poly, num_mono): dist = np.linalg.norm(r_ij, axis=1) rh.append(1. / np.mean(1. / dist)) rh = np.array(rh) - return np.mean(rh), np.std(rh) + return (np.mean(rh), 0. if num_mono == 1 else np.std(rh, ddof=0)) # python version of the espresso core function, # does not check mirror distances @@ -123,17 +131,20 @@ def test_observables_no_pbc(self): core_re = self.system.analysis.calc_re(chain_start=0, number_of_chains=num_poly, chain_length=num_mono) - np.testing.assert_allclose(core_re, self.calc_re(num_poly, num_mono)) + np.testing.assert_allclose(core_re, self.calc_re(num_poly, num_mono), + atol=1e-8) # compare calc_rg() core_rg = self.system.analysis.calc_rg(chain_start=0, number_of_chains=num_poly, chain_length=num_mono) - np.testing.assert_allclose(core_rg, self.calc_rg(num_poly, num_mono)) + np.testing.assert_allclose(core_rg, self.calc_rg(num_poly, num_mono), + atol=1e-8) # compare calc_rh() core_rh = self.system.analysis.calc_rh(chain_start=0, number_of_chains=num_poly, chain_length=num_mono) - np.testing.assert_allclose(core_rh, self.calc_rh(num_poly, num_mono)) + np.testing.assert_allclose(core_rh, self.calc_rh(num_poly, num_mono), + atol=1e-8) def test_observables_lebc(self): lin_protocol = espressomd.lees_edwards.LinearShear( @@ -163,7 +174,8 @@ def test_observables_lebc(self): core_re = self.system.analysis.calc_re(chain_start=0, number_of_chains=num_poly, chain_length=num_mono) - np.testing.assert_allclose(core_re, self.calc_re(num_poly, num_mono)) + np.testing.assert_allclose(core_re, self.calc_re(num_poly, num_mono), + atol=1e-8) # compare calc_rg() core_rg = self.system.analysis.calc_rg(chain_start=0, number_of_chains=num_poly, @@ -174,7 +186,8 @@ def test_observables_lebc(self): core_rh = self.system.analysis.calc_rh(chain_start=0, number_of_chains=num_poly, chain_length=num_mono) - np.testing.assert_allclose(core_rh, self.calc_rh(num_poly, num_mono)) + np.testing.assert_allclose(core_rh, self.calc_rh(num_poly, num_mono), + atol=1e-8) def test_exceptions(self): err_msg = """particle with id 10 does not exist @@ -185,9 +198,17 @@ def test_exceptions(self): self.insert_polymers(num_poly, num_mono) analysis = self.system.analysis for method in (analysis.calc_re, analysis.calc_rg, analysis.calc_rh): - with self.assertRaisesRegex(ValueError, err_msg): + with self.assertRaisesRegex(RuntimeError, err_msg): method(chain_start=0, number_of_chains=num_poly, chain_length=2 * num_mono) + with self.assertRaisesRegex(ValueError, "needs at least 1 bead per chain"): + method(chain_start=0, number_of_chains=1, chain_length=0) + with self.assertRaisesRegex(ValueError, "needs at least 1 bead per chain"): + method(chain_start=0, number_of_chains=1, chain_length=-1) + with self.assertRaisesRegex(ValueError, "needs at least 1 chain"): + method(chain_start=0, number_of_chains=0, chain_length=1) + with self.assertRaisesRegex(ValueError, "needs at least 1 chain"): + method(chain_start=0, number_of_chains=-1, chain_length=1) if __name__ == "__main__": diff --git a/testsuite/python/integrator_npt_stats.py b/testsuite/python/integrator_npt_stats.py index 58b209f6400..11ecb4db30c 100644 --- a/testsuite/python/integrator_npt_stats.py +++ b/testsuite/python/integrator_npt_stats.py @@ -70,7 +70,7 @@ def test_compressibility(self): compressibility = np.var(Vs) / np.average(Vs) self.assertAlmostEqual(avp, p_ext, delta=0.02) - self.assertAlmostEqual(compressibility, 0.32, delta=0.025) + self.assertAlmostEqual(compressibility, 0.32, delta=0.05) if __name__ == "__main__":