From 12b1c3a77060a0d402984ce88b230744fa31c79d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 1 Sep 2022 13:37:40 +0200 Subject: [PATCH] core: Fix broken FFT check Properly check the magnitude of complex residuals. Add a parameter to disable the check for debugging purposes. --- src/core/electrostatics/p3m.cpp | 8 +++++--- src/core/electrostatics/p3m.hpp | 3 ++- src/core/p3m/fft.cpp | 2 +- .../unit_tests/EspressoSystemStandAlone_test.cpp | 2 +- src/python/espressomd/electrostatics.py | 7 +++++++ src/script_interface/electrostatics/CoulombP3M.hpp | 6 ++++-- .../electrostatics/CoulombP3MGPU.hpp | 6 ++++-- testsuite/python/coulomb_interface.py | 6 ++++-- testsuite/python/icc_interface.py | 14 ++++++++++---- testsuite/python/save_checkpoint.py | 1 + testsuite/python/test_checkpoint.py | 2 ++ 11 files changed, 41 insertions(+), 16 deletions(-) diff --git a/src/core/electrostatics/p3m.cpp b/src/core/electrostatics/p3m.cpp index d0d059157f1..b06300fed9e 100644 --- a/src/core/electrostatics/p3m.cpp +++ b/src/core/electrostatics/p3m.cpp @@ -270,9 +270,11 @@ void CoulombP3M::init() { } CoulombP3M::CoulombP3M(P3MParameters &¶meters, double prefactor, - int tune_timings, bool tune_verbose) + int tune_timings, bool tune_verbose, + bool check_complex_residuals) : p3m{std::move(parameters)}, tune_timings{tune_timings}, - tune_verbose{tune_verbose} { + tune_verbose{tune_verbose}, check_complex_residuals{ + check_complex_residuals} { m_is_tuned = !p3m.params.tuning; p3m.params.tuning = false; @@ -490,7 +492,7 @@ double CoulombP3M::long_range_kernel(bool force_flag, bool energy_flag, } /* Back FFT force component mesh */ - auto const check_complex = !p3m.params.tuning; + auto const check_complex = !p3m.params.tuning and check_complex_residuals; for (int d = 0; d < 3; d++) { fft_perform_back(p3m.E_mesh[d].data(), check_complex, p3m.fft, comm_cart); } diff --git a/src/core/electrostatics/p3m.hpp b/src/core/electrostatics/p3m.hpp index 4f682bf7fa8..a8d02b055b4 100644 --- a/src/core/electrostatics/p3m.hpp +++ b/src/core/electrostatics/p3m.hpp @@ -89,13 +89,14 @@ struct CoulombP3M : public Coulomb::Actor { int tune_timings; bool tune_verbose; + bool check_complex_residuals; private: bool m_is_tuned; public: CoulombP3M(P3MParameters &¶meters, double prefactor, int tune_timings, - bool tune_verbose); + bool tune_verbose, bool check_complex_residuals); bool is_tuned() const { return m_is_tuned; } diff --git a/src/core/p3m/fft.cpp b/src/core/p3m/fft.cpp index 55ff03c66d3..2c1d193155c 100644 --- a/src/core/p3m/fft.cpp +++ b/src/core/p3m/fft.cpp @@ -738,7 +738,7 @@ void fft_perform_back(double *data, bool check_complex, fft_data_struct &fft, for (int i = 0; i < fft.plan[1].new_size; i++) { fft.data_buf[i] = data[2 * i]; /* real value */ // Vincent: - if (check_complex && (data[2 * i + 1] > 1e-5)) { + if (check_complex and std::abs(data[2 * i + 1]) > 1e-5) { printf("Complex value is not zero (i=%d,data=%g)!!!\n", i, data[2 * i + 1]); if (i > 100) diff --git a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp index 94cbc51a24d..ddf19a283ed 100644 --- a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp +++ b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp @@ -115,7 +115,7 @@ static void mpi_set_tuned_p3m_local(double prefactor) { 0.654, 1e-3}; auto solver = - std::make_shared(std::move(p3m), prefactor, 1, false); + std::make_shared(std::move(p3m), prefactor, 1, false, true); ::Coulomb::add_actor(solver); } diff --git a/src/python/espressomd/electrostatics.py b/src/python/espressomd/electrostatics.py index 6b91ec6d355..60a12cc4483 100644 --- a/src/python/espressomd/electrostatics.py +++ b/src/python/espressomd/electrostatics.py @@ -155,6 +155,7 @@ def default_params(self): "mesh_off": [-1., -1., -1.], "prefactor": 0., "check_neutrality": True, + "check_complex_residuals": True, "tune": True, "timings": 10, "verbose": True} @@ -231,6 +232,9 @@ class P3M(_P3MBase): check_neutrality : :obj:`bool`, optional Raise a warning if the system is not electrically neutral when set to ``True`` (default). + check_complex_residuals: :obj:`bool`, optional + Raise a warning if the backward Fourier transform has non-zero + complex residuals when set to ``True`` (default). """ _so_name = "Coulomb::CoulombP3M" @@ -281,6 +285,9 @@ class P3MGPU(_P3MBase): check_neutrality : :obj:`bool`, optional Raise a warning if the system is not electrically neutral when set to ``True`` (default). + check_complex_residuals: :obj:`bool`, optional + Raise a warning if the backward Fourier transform has non-zero + complex residuals when set to ``True`` (default). """ _so_name = "Coulomb::CoulombP3MGPU" diff --git a/src/script_interface/electrostatics/CoulombP3M.hpp b/src/script_interface/electrostatics/CoulombP3M.hpp index 147ee7d0af7..0f167cfa442 100644 --- a/src/script_interface/electrostatics/CoulombP3M.hpp +++ b/src/script_interface/electrostatics/CoulombP3M.hpp @@ -69,6 +69,8 @@ class CoulombP3M : public Actor { {"timings", AutoParameter::read_only, [this]() { return actor()->tune_timings; }}, {"tune", AutoParameter::read_only, [this]() { return m_tune; }}, + {"check_complex_residuals", AutoParameter::read_only, + [this]() { return actor()->check_complex_residuals; }}, }); } @@ -85,8 +87,8 @@ class CoulombP3M : public Actor { get_value(params, "accuracy")}; m_actor = std::make_shared( std::move(p3m), get_value(params, "prefactor"), - get_value(params, "timings"), - get_value(params, "verbose")); + get_value(params, "timings"), get_value(params, "verbose"), + get_value(params, "check_complex_residuals")); }); set_charge_neutrality_tolerance(params); } diff --git a/src/script_interface/electrostatics/CoulombP3MGPU.hpp b/src/script_interface/electrostatics/CoulombP3MGPU.hpp index d212a5e671f..3fe11a9de85 100644 --- a/src/script_interface/electrostatics/CoulombP3MGPU.hpp +++ b/src/script_interface/electrostatics/CoulombP3MGPU.hpp @@ -70,6 +70,8 @@ class CoulombP3MGPU : public Actor { {"timings", AutoParameter::read_only, [this]() { return actor()->tune_timings; }}, {"tune", AutoParameter::read_only, [this]() { return m_tune; }}, + {"check_complex_residuals", AutoParameter::read_only, + [this]() { return actor()->check_complex_residuals; }}, }); } @@ -86,8 +88,8 @@ class CoulombP3MGPU : public Actor { get_value(params, "accuracy")}; m_actor = std::make_shared( std::move(p3m), get_value(params, "prefactor"), - get_value(params, "timings"), - get_value(params, "verbose")); + get_value(params, "timings"), get_value(params, "verbose"), + get_value(params, "check_complex_residuals")); }); m_actor->request_gpu(); set_charge_neutrality_tolerance(params); diff --git a/testsuite/python/coulomb_interface.py b/testsuite/python/coulomb_interface.py index eff174e6e08..a5a1dc6b828 100644 --- a/testsuite/python/coulomb_interface.py +++ b/testsuite/python/coulomb_interface.py @@ -59,7 +59,8 @@ def tearDown(self): system, espressomd.electrostatics.P3M, dict(prefactor=2., epsilon=0., mesh_off=[0.6, 0.7, 0.8], r_cut=1.5, cao=2, mesh=[8, 10, 8], alpha=12., accuracy=0.01, tune=False, - check_neutrality=True, charge_neutrality_tolerance=7e-12)) + check_neutrality=True, charge_neutrality_tolerance=7e-12, + check_complex_residuals=False)) test_p3m_cpu_non_metallic = tests_common.generate_test_for_actor_class( system, espressomd.electrostatics.P3M, dict(prefactor=2., epsilon=3., mesh_off=[0.6, 0.7, 0.8], r_cut=1.5, @@ -78,7 +79,8 @@ def tearDown(self): system, espressomd.electrostatics.P3MGPU, dict(prefactor=2., epsilon=0., mesh_off=[0.6, 0.7, 0.8], r_cut=1.5, cao=2, mesh=[8, 10, 8], alpha=12., accuracy=0.01, tune=False, - check_neutrality=True, charge_neutrality_tolerance=7e-12)) + check_neutrality=True, charge_neutrality_tolerance=7e-12, + check_complex_residuals=False)) test_p3m_gpu_non_metallic = tests_common.generate_test_for_actor_class( system, espressomd.electrostatics.P3MGPU, dict(prefactor=2., epsilon=3., mesh_off=[0.6, 0.7, 0.8], r_cut=1.5, diff --git a/testsuite/python/icc_interface.py b/testsuite/python/icc_interface.py index 2b968d2432b..6eadca2d0e7 100644 --- a/testsuite/python/icc_interface.py +++ b/testsuite/python/icc_interface.py @@ -121,15 +121,21 @@ def test_exceptions_small_r_cut(self): @utx.skipIfMissingFeatures(["P3M"]) def test_exceptions_large_r_cut(self): - icc, (_, p) = self.setup_icc_particles_and_solver(max_iterations=1) - p3m = espressomd.electrostatics.P3M(**self.valid_p3m_parameters()) + icc, (_, p) = self.setup_icc_particles_and_solver( + max_iterations=1, convergence=10.) + p3m = espressomd.electrostatics.P3M( + check_complex_residuals=False, **self.valid_p3m_parameters()) self.system.actors.add(p3m) self.system.actors.add(icc) - with self.assertRaisesRegex(Exception, f"Particle with id {p.id} has a charge .+ that is too large for the ICC algorithm"): - p.q = 1e9 + p.q = 1e8 self.system.integrator.run(0) + + self.system.actors.remove(icc) + self.system.part.clear() + icc, (_, p) = self.setup_icc_particles_and_solver(max_iterations=1) + self.system.actors.add(icc) with self.assertRaisesRegex(Exception, "ICC failed to converge in the given number of maximal steps"): p.q = 0. self.system.integrator.run(0) diff --git a/testsuite/python/save_checkpoint.py b/testsuite/python/save_checkpoint.py index 09ebda6f002..55ab1b37b82 100644 --- a/testsuite/python/save_checkpoint.py +++ b/testsuite/python/save_checkpoint.py @@ -125,6 +125,7 @@ cao=1, alpha=1.0, r_cut=1.0, + check_complex_residuals=False, timings=15, tune=False) if 'ELC' in modes: diff --git a/testsuite/python/test_checkpoint.py b/testsuite/python/test_checkpoint.py index 17363085c3f..98f9959f68f 100644 --- a/testsuite/python/test_checkpoint.py +++ b/testsuite/python/test_checkpoint.py @@ -504,6 +504,7 @@ def test_p3m(self): reference = {'prefactor': 1.0, 'accuracy': 0.1, 'mesh': 3 * [10], 'cao': 1, 'alpha': 1.0, 'r_cut': 1.0, 'tune': False, 'timings': 15, 'check_neutrality': True, + 'check_complex_residuals': False, 'charge_neutrality_tolerance': 1e-12} for key in reference: self.assertIn(key, state) @@ -519,6 +520,7 @@ def test_elc(self): p3m_reference = {'prefactor': 1.0, 'accuracy': 0.1, 'mesh': 3 * [10], 'cao': 1, 'alpha': 1.0, 'r_cut': 1.0, 'tune': False, 'timings': 15, 'check_neutrality': True, + 'check_complex_residuals': False, 'charge_neutrality_tolerance': 7e-12} elc_reference = {'gap_size': 6.0, 'maxPWerror': 0.1, 'delta_mid_top': 0.9, 'delta_mid_bot': 0.1,