Add support for ARM architectures (espressomd#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

@@ -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}"

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<double>::infinity();
+          max_rel_diff = std::numeric_limits<double>::max();
           break;
         }
       }

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<double>::epsilon();
+  auto constexpr tol = 8. * 100. * std::numeric_limits<double>::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);
       }
     }
   }

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<double>::epsilon();
+  auto constexpr tol = 8. * 100. * std::numeric_limits<double>::epsilon();
 
   // check acceptance rate
   ReactionAlgorithmTest r_algo(42, 1., 0., {});

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<double>::epsilon();
+  auto constexpr tol = 8. * 100. * std::numeric_limits<double>::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);
         }
       }
     }

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<double>::epsilon();
+  auto constexpr tol = 8. * 100. * std::numeric_limits<double>::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<int, int>{{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);
       }
     }
   }

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<double>::epsilon();
+  auto constexpr tol = 8. * 100. * std::numeric_limits<double>::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);
     }
   }
 }

src/core/statistics_chain.cpp

@@ -34,12 +34,11 @@
 #include <cmath>
 #include <stdexcept>
 
-std::array<double, 4> calc_re(int chain_start, int chain_n_chains,
-                              int chain_length) {
+std::array<double, 4> calc_re(int chain_start, int n_chains, int chain_length) {
   double dist = 0.0, dist2 = 0.0, dist4 = 0.0;
   std::array<double, 4> 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<double, 4> calc_re(int chain_start, int chain_n_chains,
     dist2 += norm2;
     dist4 += norm2 * norm2;
   }
-  auto const tmp = static_cast<double>(chain_n_chains);
+  auto const tmp = static_cast<double>(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<double, 4> calc_rg(int chain_start, int chain_n_chains,
-                              int chain_length) {
+std::array<double, 4> 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<double, 4> 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<double, 4> calc_rg(int chain_start, int chain_n_chains,
     r_G2 += tmp;
     r_G4 += tmp * tmp;
   }
-  auto const tmp = static_cast<double>(chain_n_chains);
+  auto const tmp = static_cast<double>(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<double, 2> calc_rh(int chain_start, int chain_n_chains,
-                              int chain_length) {
+std::array<double, 2> calc_rh(int chain_start, int n_chains, int chain_length) {
   double r_H = 0.0, r_H2 = 0.0;
   std::array<double, 2> rh;
 
   auto const chain_l = static_cast<double>(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<double, 2> calc_rh(int chain_start, int chain_n_chains,
     r_H += tmp;
     r_H2 += tmp * tmp;
   }
-  auto const tmp = static_cast<double>(chain_n_chains);
+  auto const tmp = static_cast<double>(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;
 }

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<double, 4> calc_re(int chain_start, int chain_n_chains,
-                              int chain_length);
+std::array<double, 4> calc_re(int chain_start, int n_chains, int chain_length);
 
 /**
  * @brief Calculate the radius of gyration.
@@ -51,11 +50,10 @@ std::array<double, 4> 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<double, 4> calc_rg(int chain_start, int chain_n_chains,
-                              int chain_length);
+std::array<double, 4> 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<double, 4> 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<double, 2> calc_rh(int chain_start, int chain_n_chains,
-                              int chain_length);
+std::array<double, 2> calc_rh(int chain_start, int n_chains, int chain_length);
 /**@}*/
 
 #endif

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<double>::epsilon();
+  auto constexpr tol = 8. * 100. * std::numeric_limits<double>::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

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<double>::epsilon();
+  auto constexpr tol = 8. * 100. * std::numeric_limits<double>::epsilon();
   boost::mpi::communicator world;
 
   auto const box_l = 8.;

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 <http://www.gnu.org/licenses/>.
  */
-#define BOOST_TEST_MODULE AutoParameter test
+
+#define BOOST_TEST_MODULE Field coupling test for fields
 #define BOOST_TEST_DYN_LINK
 #include <boost/test/unit_test.hpp>
 

src/core/unit_tests/field_coupling_fields_test.cpp

@@ -17,7 +17,7 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#define BOOST_TEST_MODULE AutoParameter test
+#define BOOST_TEST_MODULE Field coupling test for fields
 #define BOOST_TEST_DYN_LINK
 #include <boost/test/unit_test.hpp>
 
@@ -37,7 +37,7 @@
 #include <limits>
 #include <type_traits>
 
-constexpr auto eps = 10 * std::numeric_limits<double>::epsilon();
+auto constexpr eps = 10. * std::numeric_limits<double>::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__
 }

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 <http://www.gnu.org/licenses/>.
  */
-#define BOOST_TEST_MODULE test
+
+#define BOOST_TEST_MODULE Field coupling test for force fields
 #define BOOST_TEST_DYN_LINK
 #include <boost/test/unit_test.hpp>
 

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<double>::epsilon();
+  auto constexpr tol = 8. * 100. * std::numeric_limits<double>::epsilon();
 
   // check only trivial cases
   for (auto const cao : {1, 2, 3, 4, 5, 6, 7}) {

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<double>::epsilon());
+                4. * std::numeric_limits<double>::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<double>::epsilon());
+                4. * std::numeric_limits<double>::epsilon());
   }
 }
 

src/core/unit_tests/specfunc_test.cpp

@@ -26,7 +26,7 @@
 #include <cmath>
 #include <limits>
 
-constexpr auto eps = 100. * std::numeric_limits<double>::epsilon();
+auto constexpr eps = 8. * 100. * std::numeric_limits<double>::epsilon();
 
 BOOST_AUTO_TEST_CASE(hurwitz_zeta_function) {
   constexpr auto max_bits = 54.0;

src/core/unit_tests/thermostats_test.cpp

@@ -41,8 +41,8 @@
 #include <tuple>
 
 // 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<double>::epsilon();
+// and by 8 to account for error accumulation in thermostat functions
+auto constexpr tol = 8. * 100. * std::numeric_limits<double>::epsilon();
 
 Particle particle_factory() {
   Particle p{};

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<double>::epsilon();
+  auto constexpr tol = 8. * 100. * std::numeric_limits<double>::epsilon();
   std::vector<double> const values{1., 2., 3., 4.};
   {
     auto const res = WeightedAverage<Testing::Identity, Testing::One>()(values);