Improve resilience to precision loss

.github/workflows/push_pull.yml

@@ -40,6 +40,9 @@ jobs:
           build_procs: 3
           check_procs: 3
           with_ccache: 'true'
+          with_walberla: 'true'
+          with_walberla_fft: 'false'
+          with_walberla_avx: 'false'
 
   debian:
     runs-on: ubuntu-latest

maintainer/CI/build_cmake.sh

@@ -132,6 +132,7 @@ set_default_value with_fftw true
 set_default_value with_gsl true
 set_default_value with_scafacos false
 set_default_value with_walberla false
+set_default_value with_walberla_fft true
 set_default_value with_walberla_avx false
 set_default_value with_stokesian_dynamics false
 set_default_value test_timeout 500
@@ -168,7 +169,9 @@ cmake_params="${cmake_params} -D ESPRESSO_BUILD_WITH_STOKESIAN_DYNAMICS=${with_s
 cmake_params="${cmake_params} -D ESPRESSO_BUILD_WITH_WALBERLA=${with_walberla}"
 
 if [ "${with_walberla}" = true ]; then
-  cmake_params="${cmake_params} -D ESPRESSO_BUILD_WITH_WALBERLA_FFT=ON"
+  if [ "${with_walberla_fft}" = true ]; then
+    cmake_params="${cmake_params} -D ESPRESSO_BUILD_WITH_WALBERLA_FFT=ON"
+  fi
   if [ "${with_walberla_avx}" = true ]; then
     cmake_params="${cmake_params} -D ESPRESSO_BUILD_WITH_WALBERLA_AVX=ON"
   fi

src/utils/include/utils/Histogram.hpp

@@ -105,13 +105,10 @@ class Histogram {
       throw std::invalid_argument("Wrong dimensions for the value");
     }
     if (check_limits(pos)) {
-      boost::array<array_index, M + 1> index;
-      for (std::size_t i = 0; i < M; ++i) {
-        index[i] = calc_bin_index(pos[i], m_limits[i].first, m_bin_sizes[i]);
-      }
-      for (array_index i = 0; i < static_cast<array_index>(N); ++i) {
-        index.back() = i;
-        m_array(index) += value[static_cast<std::size_t>(i)];
+      auto index = calc_bin_index(pos);
+      for (std::size_t i = 0; i < N; ++i) {
+        index.back() = static_cast<array_index>(i);
+        m_array(index) += value[i];
         m_count(index)++;
       }
     }
@@ -121,20 +118,29 @@ class Histogram {
   virtual void normalize() {
     auto const bin_volume = std::accumulate(
         m_bin_sizes.begin(), m_bin_sizes.end(), U{1}, std::multiplies<U>());
-    std::transform(
-        m_array.data(), m_array.data() + m_array.num_elements(), m_array.data(),
+    std::ranges::transform(
+        std::span(m_array.data(), m_array.num_elements()), m_array.data(),
         [bin_volume](T v) { return static_cast<T>(v / bin_volume); });
   }
 
 private:
   /**
    * \brief Calculate the bin index.
-   * \param value  Position on that dimension.
-   * \param offset Bin offset on that dimension.
-   * \param size   Bin size on that dimension.
+   * \param pos  Position.
    */
-  array_index calc_bin_index(double value, double offset, double size) const {
-    return static_cast<array_index>(std::floor((value - offset) / size));
+  auto calc_bin_index(std::span<const U> const &pos) const {
+    boost::array<array_index, M + 1> index;
+    for (std::size_t i = 0; i < M; ++i) {
+      auto const offset = m_limits[i].first;
+      auto const size = m_bin_sizes[i];
+      auto const n_bins = static_cast<long>(m_n_bins[i]);
+      auto const bin = static_cast<long>(std::floor((pos[i] - offset) / size));
+      // handle edge cases when the position is exactly between two bins:
+      // due to precision loss in the offset subtraction, the bin index might
+      // be off by one, so we fold it here back inside the valid range
+      index[i] = static_cast<array_index>(std::clamp(bin, 0l, n_bins - 1l));
+    }
+    return index;
   }
 
   /**
@@ -153,19 +159,19 @@ class Histogram {
    * \brief Check if the position lies within the histogram limits.
    * \param pos     Position to check.
    */
-  bool check_limits(std::span<const U> pos) const {
+  bool check_limits(std::span<const U> const &pos) const {
     assert(pos.size() == M);
-    bool within_range = true;
-    for (std::size_t i = 0; i < M; ++i) {
-      if (pos[i] < m_limits[i].first or pos[i] >= m_limits[i].second)
-        within_range = false;
-    }
-    return within_range;
+    auto it_limits = m_limits.begin();
+    return std::ranges::all_of(pos, [&it_limits](U const value) {
+      auto const [lower, upper] = *it_limits;
+      ++it_limits;
+      return value >= lower and value <= upper;
+    });
   }
 
   std::array<std::size_t, M + 1> m_array_dim() const {
     std::array<std::size_t, M + 1> dimensions;
-    std::copy(m_n_bins.begin(), m_n_bins.end(), dimensions.begin());
+    std::ranges::copy(m_n_bins, dimensions.begin());
     dimensions.back() = N;
     return dimensions;
   }
@@ -193,14 +199,15 @@ class Histogram {
  */
 template <typename T, std::size_t N, std::size_t M = 3, typename U = double>
 class CylindricalHistogram : public Histogram<T, N, M, U> {
-  using Histogram<T, N, M, U>::m_n_bins;
-  using Histogram<T, N, M, U>::m_limits;
-  using Histogram<T, N, M, U>::m_bin_sizes;
-  using Histogram<T, N, M, U>::m_array;
-  using typename Histogram<T, N, M, U>::array_index;
+  using Base = Histogram<T, N, M, U>;
+  using Base::m_array;
+  using Base::m_bin_sizes;
+  using Base::m_limits;
+  using Base::m_n_bins;
+  using typename Base::array_index;
 
 public:
-  using Histogram<T, N, M, U>::Histogram;
+  using Base::Histogram;
 
   void normalize() override {
     auto const min_r = m_limits[0].first;
@@ -214,8 +221,8 @@ class CylindricalHistogram : public Histogram<T, N, M, U> {
       auto const bin_volume = (r_right * r_right - r_left * r_left) *
                               z_bin_size * phi_bin_size / U(2);
       auto *begin = m_array[i].origin();
-      std::transform(
-          begin, begin + m_array[i].num_elements(), begin,
+      std::ranges::transform(
+          std::span(begin, m_array[i].num_elements()), begin,
           [bin_volume](T v) { return static_cast<T>(v / bin_volume); });
     }
   }

testsuite/python/elc_vs_analytic.py

@@ -114,7 +114,7 @@ class TestCPU(Test, ut.TestCase):
 class TestGPU(Test, ut.TestCase):
 
     p3m_class = espressomd.electrostatics.P3MGPU
-    rtol = 4e-6
+    rtol = 5e-6
 
 
 if __name__ == "__main__":

testsuite/python/lb.py

@@ -852,8 +852,8 @@ class LBTestWalberlaSinglePrecisionCPU(LBTest, ut.TestCase):
     lb_class = espressomd.lb.LBFluidWalberla
     lb_lattice_class = espressomd.lb.LatticeWalberla
     lb_params = {"single_precision": True}
-    atol = 1e-7
-    rtol = 5e-5
+    atol = 5e-6
+    rtol = 2e-4
 
 
 @utx.skipIfMissingGPU()
@@ -872,7 +872,7 @@ class LBTestWalberlaSinglePrecisionGPU(LBTest, ut.TestCase):
     lb_class = espressomd.lb.LBFluidWalberlaGPU
     lb_lattice_class = espressomd.lb.LatticeWalberla
     lb_params = {"single_precision": True}
-    atol = 1e-6
+    atol = 5e-6
     rtol = 2e-4
 
 
@@ -894,7 +894,7 @@ class LBTestWalberlaSinglePrecisionBlocksCPU(LBTest, ut.TestCase):
     blocks_per_mpi_rank = [2, 2, 2]
     lb_params = {"single_precision": True,
                  "blocks_per_mpi_rank": blocks_per_mpi_rank}
-    atol = 1e-6
+    atol = 5e-6
     rtol = 2e-4
 
 

testsuite/python/lb_electrohydrodynamics.py

@@ -71,7 +71,7 @@ def test(self):
 
         system.integrator.run(steps=500)
 
-        np.testing.assert_allclose(v_term, np.copy(p.v), atol=5e-5)
+        np.testing.assert_allclose(v_term, np.copy(p.v), atol=6e-5)
 
 
 if __name__ == "__main__":

testsuite/python/lb_lees_edwards_particle_coupling.py

@@ -328,7 +328,7 @@ def test_momentum_conservation(self):
                 np.sum(lbf[:, :, :].last_applied_force, axis=(0, 1, 2))), atol=1E-9)
             current_mom = np.copy(system.analysis.linear_momentum())
             np.testing.assert_allclose(
-                initial_mom[1:], current_mom[1:], atol=2.75E-7)
+                initial_mom[1:], current_mom[1:], atol=1E-6)
 
 
 if __name__ == '__main__':

testsuite/python/lb_planar_couette.py

@@ -107,7 +107,7 @@ def check_profile(self, u_getter, **kwargs):
             u_ref = analytical(pos, system.time - 1. + 1., lbf.kinematic_viscosity,
                                shear_velocity, h, k_max)
             u_lbf = np.copy(u_getter(lbf).reshape([-1]))
-            np.testing.assert_allclose(u_lbf, u_ref, atol=1e-4, rtol=0.)
+            np.testing.assert_allclose(u_lbf, u_ref, atol=2e-4, rtol=0.)
 
     def test_profile_xy(self):
         if "blocks_per_mpi_rank" in self.lb_params:

testsuite/python/oif_volume_conservation.py

@@ -90,7 +90,7 @@ def check_relaxation(self, **kwargs):
             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)
+        self.assertAlmostEqual(lam, 325e-4, delta=5e-4)
         self.system.thermostat.turn_off()
         self.system.part.clear()
 

testsuite/python/p3m_fft.py

@@ -103,7 +103,7 @@ def test_fft_plans(self):
                 ref_energy = -75.871906
                 p3m_energy = self.system.analysis.energy()['coulomb']
                 self.system.electrostatics.clear()
-                np.testing.assert_allclose(p3m_energy, ref_energy, rtol=1e-4)
+                np.testing.assert_allclose(p3m_energy, ref_energy, rtol=1e-3)
 
     @utx.skipIfMissingFeatures("DP3M")
     @ut.skipIf(n_nodes < 2 or n_nodes >= 8, "only runs for 2 <= n_nodes <= 7")

testsuite/python/p3m_madelung.py

@@ -215,7 +215,7 @@ def test_infinite_ionic_wire(self):
 
         def check():
             energy, p_scalar, p_tensor = self.get_normalized_obs_per_ion()
-            np.testing.assert_allclose(energy, ref_energy, atol=0., rtol=5e-7)
+            np.testing.assert_allclose(energy, ref_energy, atol=0., rtol=7e-7)
             np.testing.assert_allclose(p_scalar, np.trace(ref_pressure) / 3.,
                                        atol=1e-12, rtol=1e-6)
             np.testing.assert_allclose(p_tensor, ref_pressure, atol=1e-12,
@@ -291,7 +291,7 @@ def test_infinite_ionic_cube(self):
 
         def check():
             energy, p_scalar, p_tensor = self.get_normalized_obs_per_ion()
-            np.testing.assert_allclose(energy, ref_energy, atol=0., rtol=1e-6)
+            np.testing.assert_allclose(energy, ref_energy, atol=0., rtol=2e-6)
             np.testing.assert_allclose(p_scalar, np.trace(ref_pressure) / 3.,
                                        atol=1e-12, rtol=5e-6)
             np.testing.assert_allclose(p_tensor, ref_pressure, atol=5e-9,