diff --git a/.github/workflows/push_pull.yml b/.github/workflows/push_pull.yml index 233254e4aa0..12a0bf74eb2 100644 --- a/.github/workflows/push_pull.yml +++ b/.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 diff --git a/maintainer/CI/build_cmake.sh b/maintainer/CI/build_cmake.sh index 0e653ebbe5f..0c8f7314e62 100755 --- a/maintainer/CI/build_cmake.sh +++ b/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 diff --git a/src/core/observables/RDF.cpp b/src/core/observables/RDF.cpp index a42fd7d8c84..1ae6a7e9859 100644 --- a/src/core/observables/RDF.cpp +++ b/src/core/observables/RDF.cpp @@ -29,6 +29,7 @@ #include #include +#include #include #include #include @@ -71,8 +72,8 @@ RDF::evaluate(boost::mpi::communicator const &comm, auto const &pos2) { auto const dist = box_geo.get_mi_vector(pos1, pos2).norm(); if (dist > min_r && dist < max_r) { - auto const ind = - static_cast(std::floor((dist - min_r) * inv_bin_width)); + auto ind = static_cast(std::floor((dist - min_r) * inv_bin_width)); + ind = std::clamp(ind, 0l, static_cast(n_r_bins) - 1l); res[ind]++; } cnt++; diff --git a/src/python/espressomd/observables.py b/src/python/espressomd/observables.py index 49dc8714899..8b6b457f470 100644 --- a/src/python/espressomd/observables.py +++ b/src/python/espressomd/observables.py @@ -75,8 +75,8 @@ class CylindricalProfileObservable(ProfileObservable): Base class for observables that work with cylinder coordinates """ - def __init__( - self, transform_params=CylindricalTransformationParameters(), **kwargs): + def __init__(self, transform_params=CylindricalTransformationParameters(), + **kwargs): # Provide default transformation parameters if not user-provided kwargs['transform_params'] = transform_params super().__init__(**kwargs) @@ -152,17 +152,17 @@ class DensityProfile(ProfileObservable): n_z_bins : :obj:`int` Number of bins in ``z`` direction. min_x : :obj:`float` - Minimum ``x`` to consider. + Minimum ``x`` to consider (inclusive). min_y : :obj:`float` - Minimum ``y`` to consider. + Minimum ``y`` to consider (inclusive). min_z : :obj:`float` - Minimum ``z`` to consider. + Minimum ``z`` to consider (inclusive). max_x : :obj:`float` - Maximum ``x`` to consider. + Maximum ``x`` to consider (exclusive). max_y : :obj:`float` - Maximum ``y`` to consider. + Maximum ``y`` to consider (exclusive). max_z : :obj:`float` - Maximum ``z`` to consider. + Maximum ``z`` to consider (exclusive). Methods ------- @@ -218,17 +218,17 @@ class FluxDensityProfile(ProfileObservable): n_z_bins : :obj:`int` Number of bins in ``z`` direction. min_x : :obj:`float` - Minimum ``x`` to consider. + Minimum ``x`` to consider (inclusive). min_y : :obj:`float` - Minimum ``y`` to consider. + Minimum ``y`` to consider (inclusive). min_z : :obj:`float` - Minimum ``z`` to consider. + Minimum ``z`` to consider (inclusive). max_x : :obj:`float` - Maximum ``x`` to consider. + Maximum ``x`` to consider (exclusive). max_y : :obj:`float` - Maximum ``y`` to consider. + Maximum ``y`` to consider (exclusive). max_z : :obj:`float` - Maximum ``z`` to consider. + Maximum ``z`` to consider (exclusive). Methods ------- @@ -261,17 +261,17 @@ class ForceDensityProfile(ProfileObservable): n_z_bins : :obj:`int` Number of bins in ``z`` direction. min_x : :obj:`float` - Minimum ``x`` to consider. + Minimum ``x`` to consider (inclusive). min_y : :obj:`float` - Minimum ``y`` to consider. + Minimum ``y`` to consider (inclusive). min_z : :obj:`float` - Minimum ``z`` to consider. + Minimum ``z`` to consider (inclusive). max_x : :obj:`float` - Maximum ``x`` to consider. + Maximum ``x`` to consider (exclusive). max_y : :obj:`float` - Maximum ``y`` to consider. + Maximum ``y`` to consider (exclusive). max_z : :obj:`float` - Maximum ``z`` to consider. + Maximum ``z`` to consider (exclusive). Methods ------- @@ -306,17 +306,17 @@ class LBVelocityProfile(ProfileObservable): n_z_bins : :obj:`int` Number of bins in ``z`` direction. min_x : :obj:`float` - Minimum ``x`` to consider. + Minimum ``x`` to consider (inclusive). min_y : :obj:`float` - Minimum ``y`` to consider. + Minimum ``y`` to consider (inclusive). min_z : :obj:`float` - Minimum ``z`` to consider. + Minimum ``z`` to consider (inclusive). max_x : :obj:`float` - Maximum ``x`` to consider. + Maximum ``x`` to consider (exclusive). max_y : :obj:`float` - Maximum ``y`` to consider. + Maximum ``y`` to consider (exclusive). max_z : :obj:`float` - Maximum ``z`` to consider. + Maximum ``z`` to consider (exclusive). sampling_delta_x : :obj:`float`, default=1.0 Spacing for the sampling grid in ``x``-direction. sampling_delta_y : :obj:`float`, default=1.0 @@ -866,17 +866,17 @@ class CylindricalDensityProfile(CylindricalProfileObservable): n_z_bins : :obj:`int`, default = 1 Number of bins in ``z`` direction. min_r : :obj:`float`, default = 0 - Minimum ``r`` to consider. + Minimum ``r`` to consider (inclusive). min_phi : :obj:`float`, default = :math:`-\\pi` - Minimum ``phi`` to consider. Must be in :math:`[-\\pi,\\pi)`. + Minimum ``phi`` to consider (inclusive). Must be in :math:`[-\\pi,\\pi)`. min_z : :obj:`float` - Minimum ``z`` to consider. + Minimum ``z`` to consider (inclusive). max_r : :obj:`float` - Maximum ``r`` to consider. + Maximum ``r`` to consider (exclusive). max_phi : :obj:`float`, default = :math:`\\pi` - Maximum ``phi`` to consider. Must be in :math:`(-\\pi,\\pi]`. + Maximum ``phi`` to consider (exclusive). Must be in :math:`(-\\pi,\\pi]`. max_z : :obj:`float` - Maximum ``z`` to consider. + Maximum ``z`` to consider (exclusive). Methods ------- @@ -909,17 +909,17 @@ class CylindricalFluxDensityProfile(CylindricalProfileObservable): n_z_bins : :obj:`int`, default = 1 Number of bins in ``z`` direction. min_r : :obj:`float`, default = 0 - Minimum ``r`` to consider. + Minimum ``r`` to consider (inclusive). min_phi : :obj:`float`, default = :math:`-\\pi` - Minimum ``phi`` to consider. Must be in :math:`[-\\pi,\\pi)`. + Minimum ``phi`` to consider (inclusive). Must be in :math:`[-\\pi,\\pi)`. min_z : :obj:`float` - Minimum ``z`` to consider. + Minimum ``z`` to consider (inclusive). max_r : :obj:`float` - Maximum ``r`` to consider. + Maximum ``r`` to consider (exclusive). max_phi : :obj:`float`, default = :math:`\\pi` - Maximum ``phi`` to consider. Must be in :math:`(-\\pi,\\pi]`. + Maximum ``phi`` to consider (exclusive). Must be in :math:`(-\\pi,\\pi]`. max_z : :obj:`float` - Maximum ``z`` to consider. + Maximum ``z`` to consider (exclusive). Methods ------- @@ -957,17 +957,17 @@ class CylindricalLBFluxDensityProfileAtParticlePositions( n_z_bins : :obj:`int`, default = 1 Number of bins in ``z`` direction. min_r : :obj:`float`, default = 0 - Minimum ``r`` to consider. + Minimum ``r`` to consider (inclusive). min_phi : :obj:`float`, default = :math:`-\\pi` - Minimum ``phi`` to consider. Must be in :math:`[-\\pi,\\pi)`. + Minimum ``phi`` to consider (inclusive). Must be in :math:`[-\\pi,\\pi)`. min_z : :obj:`float` - Minimum ``z`` to consider. + Minimum ``z`` to consider (inclusive). max_r : :obj:`float` - Maximum ``r`` to consider. + Maximum ``r`` to consider (exclusive). max_phi : :obj:`float`, default = :math:`\\pi` - Maximum ``phi`` to consider. Must be in :math:`(-\\pi,\\pi]`. + Maximum ``phi`` to consider (exclusive). Must be in :math:`(-\\pi,\\pi]`. max_z : :obj:`float` - Maximum ``z`` to consider. + Maximum ``z`` to consider (exclusive). Methods ------- @@ -1005,17 +1005,17 @@ class CylindricalLBVelocityProfileAtParticlePositions( n_z_bins : :obj:`int`, default = 1 Number of bins in ``z`` direction. min_r : :obj:`float`, default = 0 - Minimum ``r`` to consider. + Minimum ``r`` to consider (inclusive). min_phi : :obj:`float`, default = :math:`-\\pi` - Minimum ``phi`` to consider. Must be in :math:`[-\\pi,\\pi)`. + Minimum ``phi`` to consider (inclusive). Must be in :math:`[-\\pi,\\pi)`. min_z : :obj:`float` - Minimum ``z`` to consider. + Minimum ``z`` to consider (inclusive). max_r : :obj:`float` - Maximum ``r`` to consider. + Maximum ``r`` to consider (exclusive). max_phi : :obj:`float`, default = :math:`\\pi` - Maximum ``phi`` to consider. Must be in :math:`(-\\pi,\\pi]`. + Maximum ``phi`` to consider (exclusive). Must be in :math:`(-\\pi,\\pi]`. max_z : :obj:`float` - Maximum ``z`` to consider. + Maximum ``z`` to consider (exclusive). Methods ------- @@ -1051,17 +1051,17 @@ class CylindricalVelocityProfile(CylindricalProfileObservable): n_z_bins : :obj:`int`, default = 1 Number of bins in ``z`` direction. min_r : :obj:`float`, default = 0 - Minimum ``r`` to consider. + Minimum ``r`` to consider (inclusive). min_phi : :obj:`float`, default = :math:`-\\pi` - Minimum ``phi`` to consider. Must be in :math:`[-\\pi,\\pi)`. + Minimum ``phi`` to consider (inclusive). Must be in :math:`[-\\pi,\\pi)`. min_z : :obj:`float` - Minimum ``z`` to consider. + Minimum ``z`` to consider (inclusive). max_r : :obj:`float` - Maximum ``r`` to consider. + Maximum ``r`` to consider (exclusive). max_phi : :obj:`float`, default = :math:`\\pi` - Maximum ``phi`` to consider. Must be in :math:`(-\\pi,\\pi]`. + Maximum ``phi`` to consider (exclusive). Must be in :math:`(-\\pi,\\pi]`. max_z : :obj:`float` - Maximum ``z`` to consider. + Maximum ``z`` to consider (exclusive). Methods ------- @@ -1099,17 +1099,17 @@ class CylindricalLBVelocityProfile(CylindricalProfileObservable): n_z_bins : :obj:`int`, default = 1 Number of bins in ``z`` direction. min_r : :obj:`float`, default = 0 - Minimum ``r`` to consider. + Minimum ``r`` to consider (inclusive). min_phi : :obj:`float`, default = :math:`-\\pi` - Minimum ``phi`` to consider. Must be in :math:`[-\\pi,\\pi)`. + Minimum ``phi`` to consider (inclusive). Must be in :math:`[-\\pi,\\pi)`. min_z : :obj:`float` - Minimum ``z`` to consider. + Minimum ``z`` to consider (inclusive). max_r : :obj:`float` - Maximum ``r`` to consider. + Maximum ``r`` to consider (exclusive). max_phi : :obj:`float`, default = :math:`\\pi` - Maximum ``phi`` to consider. Must be in :math:`(-\\pi,\\pi]`. + Maximum ``phi`` to consider (exclusive). Must be in :math:`(-\\pi,\\pi]`. max_z : :obj:`float` - Maximum ``z`` to consider. + Maximum ``z`` to consider (exclusive). sampling_density : :obj:`float` Samples per unit volume for the LB velocity interpolation. @@ -1145,9 +1145,9 @@ class RDF(Observable): n_r_bins : :obj:`int` Number of bins in radial direction. min_r : :obj:`float` - Minimum ``r`` to consider. + Minimum ``r`` to consider (exclusive). max_r : :obj:`float` - Maximum ``r`` to consider. + Maximum ``r`` to consider (exclusive). Methods ------- diff --git a/src/utils/include/utils/Histogram.hpp b/src/utils/include/utils/Histogram.hpp index d36421a2fb9..1623dc0ace3 100644 --- a/src/utils/include/utils/Histogram.hpp +++ b/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 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(N); ++i) { - index.back() = i; - m_array(index) += value[static_cast(i)]; + auto index = calc_bin_index(pos); + for (std::size_t i = 0; i < N; ++i) { + index.back() = static_cast(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()); - 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(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(std::floor((value - offset) / size)); + auto calc_bin_index(std::span const &pos) const { + boost::array 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(m_n_bins[i]); + auto const bin = static_cast(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(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 pos) const { + bool check_limits(std::span 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 m_array_dim() const { std::array 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 class CylindricalHistogram : public Histogram { - using Histogram::m_n_bins; - using Histogram::m_limits; - using Histogram::m_bin_sizes; - using Histogram::m_array; - using typename Histogram::array_index; + using Base = Histogram; + 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::Histogram; + using Base::Histogram; void normalize() override { auto const min_r = m_limits[0].first; @@ -214,8 +221,8 @@ class CylindricalHistogram : public Histogram { 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(v / bin_volume); }); } } diff --git a/testsuite/python/elc_vs_analytic.py b/testsuite/python/elc_vs_analytic.py index d7cb4edb722..d702821a9af 100644 --- a/testsuite/python/elc_vs_analytic.py +++ b/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__": diff --git a/testsuite/python/lb.py b/testsuite/python/lb.py index ed94748f926..444b61b6529 100644 --- a/testsuite/python/lb.py +++ b/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 diff --git a/testsuite/python/lb_electrohydrodynamics.py b/testsuite/python/lb_electrohydrodynamics.py index e8d7808f3c2..37a00f6c657 100644 --- a/testsuite/python/lb_electrohydrodynamics.py +++ b/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__": diff --git a/testsuite/python/lb_lees_edwards_particle_coupling.py b/testsuite/python/lb_lees_edwards_particle_coupling.py index ed8fbcc3328..80ab75aa584 100644 --- a/testsuite/python/lb_lees_edwards_particle_coupling.py +++ b/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__': diff --git a/testsuite/python/lb_planar_couette.py b/testsuite/python/lb_planar_couette.py index d4efc40964f..9990838a25e 100644 --- a/testsuite/python/lb_planar_couette.py +++ b/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: diff --git a/testsuite/python/oif_volume_conservation.py b/testsuite/python/oif_volume_conservation.py index ca1b0673e7a..a8c4ecbbedb 100644 --- a/testsuite/python/oif_volume_conservation.py +++ b/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() diff --git a/testsuite/python/p3m_fft.py b/testsuite/python/p3m_fft.py index 8bfb7fca51a..358f71d2dec 100644 --- a/testsuite/python/p3m_fft.py +++ b/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") diff --git a/testsuite/python/p3m_madelung.py b/testsuite/python/p3m_madelung.py index 95f263910a6..8975c0617a8 100644 --- a/testsuite/python/p3m_madelung.py +++ b/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,