Add new P3M tuning options

src/core/electrostatics/p3m.cpp

@@ -598,14 +598,17 @@ class CoulombTuningAlgorithm : public TuningAlgorithm {
   double m_mesh_density_min = -1., m_mesh_density_max = -1.;
   // indicates if mesh should be tuned
   bool m_tune_mesh = false;
+  std::pair<std::optional<int>, std::optional<int>> m_tune_limits;
 
 protected:
   P3MParameters &get_params() override { return p3m.params; }
 
 public:
   CoulombTuningAlgorithm(System::System &system, auto &input_p3m,
-                         double prefactor, int timings)
-      : TuningAlgorithm(system, prefactor, timings), p3m{input_p3m} {}
+                         double prefactor, int timings,
+                         decltype(m_tune_limits) tune_limits)
+      : TuningAlgorithm(system, prefactor, timings), p3m{input_p3m},
+        m_tune_limits{std::move(tune_limits)} {}
 
   void on_solver_change() const override { m_system.on_coulomb_change(); }
 
@@ -722,6 +725,16 @@ class CoulombTuningAlgorithm : public TuningAlgorithm {
           max_npart_per_dim, std::cbrt(static_cast<double>(p3m.sum_qpart)));
       m_mesh_density_min = min_npart_per_dim / normalized_box_dim;
       m_mesh_density_max = max_npart_per_dim / normalized_box_dim;
+      if (m_tune_limits.first or m_tune_limits.second) {
+        auto const &box_l = box_geo.length();
+        auto const dim = std::max({box_l[0], box_l[1], box_l[2]});
+        if (m_tune_limits.first) {
+          m_mesh_density_min = static_cast<double>(*m_tune_limits.first) / dim;
+        }
+        if (m_tune_limits.second) {
+          m_mesh_density_max = static_cast<double>(*m_tune_limits.second) / dim;
+        }
+      }
       m_tune_mesh = true;
     } else {
       m_mesh_density_min = m_mesh_density_max = mesh_density;
@@ -805,7 +818,7 @@ void CoulombP3MImpl<FloatType, Architecture>::tune() {
     }
     try {
       CoulombTuningAlgorithm<FloatType, Architecture> parameters(
-          system, p3m, prefactor, tune_timings);
+          system, p3m, prefactor, tune_timings, tune_limits);
       parameters.setup_logger(tune_verbose);
       // parameter ranges
       parameters.determine_mesh_limits();

src/core/electrostatics/p3m.impl.hpp

@@ -36,6 +36,7 @@
 #include <utils/Vector.hpp>
 
 #include <memory>
+#include <optional>
 #include <stdexcept>
 #include <type_traits>
 #include <utility>
@@ -69,6 +70,7 @@ struct CoulombP3MImpl : public CoulombP3M {
   /** @brief Coulomb P3M meshes and FFT algorithm. */
   std::unique_ptr<p3m_data_struct_coulomb<FloatType>> p3m_impl;
   int tune_timings;
+  std::pair<std::optional<int>, std::optional<int>> tune_limits;
   bool tune_verbose;
   bool check_complex_residuals;
   bool m_is_tuned;
@@ -77,10 +79,10 @@ struct CoulombP3MImpl : public CoulombP3M {
   CoulombP3MImpl(
       std::unique_ptr<p3m_data_struct_coulomb<FloatType>> &&p3m_handle,
       double prefactor, int tune_timings, bool tune_verbose,
-      bool check_complex_residuals)
+      decltype(tune_limits) tune_limits, bool check_complex_residuals)
       : CoulombP3M(p3m_handle->params), p3m{*p3m_handle},
         p3m_impl{std::move(p3m_handle)}, tune_timings{tune_timings},
-        tune_verbose{tune_verbose},
+        tune_limits{std::move(tune_limits)}, tune_verbose{tune_verbose},
         check_complex_residuals{check_complex_residuals} {
 
     if (tune_timings <= 0) {

src/core/magnetostatics/dp3m.cpp

@@ -74,6 +74,7 @@
 #include <sstream>
 #include <stdexcept>
 #include <tuple>
+#include <utility>
 #include <vector>
 
 #ifdef FFTW3_H
@@ -534,11 +535,14 @@ template <typename FloatType, Arch Architecture>
 class DipolarTuningAlgorithm : public TuningAlgorithm {
   p3m_data_struct_dipoles<FloatType> &dp3m;
   int m_mesh_max = -1, m_mesh_min = -1;
+  std::pair<std::optional<int>, std::optional<int>> m_tune_limits;
 
 public:
   DipolarTuningAlgorithm(System::System &system, decltype(dp3m) &input_dp3m,
-                         double prefactor, int timings)
-      : TuningAlgorithm(system, prefactor, timings), dp3m{input_dp3m} {}
+                         double prefactor, int timings,
+                         decltype(m_tune_limits) tune_limits)
+      : TuningAlgorithm(system, prefactor, timings), dp3m{input_dp3m},
+        m_tune_limits{std::move(tune_limits)} {}
 
   P3MParameters &get_params() override { return dp3m.params; }
 
@@ -603,6 +607,12 @@ class DipolarTuningAlgorithm : public TuningAlgorithm {
       m_mesh_min = static_cast<int>(std::round(std::pow(2., std::floor(expo))));
       /* avoid using more than 1 GB of FFT arrays */
       m_mesh_max = 128;
+      if (m_tune_limits.first) {
+        m_mesh_min = *m_tune_limits.first;
+      }
+      if (m_tune_limits.second) {
+        m_mesh_max = *m_tune_limits.second;
+      }
     } else {
       m_mesh_min = m_mesh_max = dp3m.params.mesh[0];
       m_logger->report_fixed_mesh(dp3m.params.mesh);
@@ -662,7 +672,7 @@ void DipolarP3MImpl<FloatType, Architecture>::tune() {
     }
     try {
       DipolarTuningAlgorithm<FloatType, Architecture> parameters(
-          system, dp3m, prefactor, tune_timings);
+          system, dp3m, prefactor, tune_timings, tune_limits);
       parameters.setup_logger(tune_verbose);
       // parameter ranges
       parameters.determine_mesh_limits();

src/core/magnetostatics/dp3m.impl.hpp

@@ -36,6 +36,7 @@
 #include <utils/Vector.hpp>
 
 #include <memory>
+#include <optional>
 #include <stdexcept>
 #include <type_traits>
 #include <utility>
@@ -72,16 +73,18 @@ struct DipolarP3MImpl : public DipolarP3M {
   /** @brief Coulomb P3M meshes and FFT algorithm. */
   std::unique_ptr<p3m_data_struct_dipoles<FloatType>> dp3m_impl;
   int tune_timings;
+  std::pair<std::optional<int>, std::optional<int>> tune_limits;
   bool tune_verbose;
   bool m_is_tuned;
 
 public:
   DipolarP3MImpl(
       std::unique_ptr<p3m_data_struct_dipoles<FloatType>> &&dp3m_handle,
-      double prefactor, int tune_timings, bool tune_verbose)
+      double prefactor, int tune_timings, bool tune_verbose,
+      decltype(tune_limits) tune_limits)
       : DipolarP3M(dp3m_handle->params), dp3m{*dp3m_handle},
         dp3m_impl{std::move(dp3m_handle)}, tune_timings{tune_timings},
-        tune_verbose{tune_verbose} {
+        tune_limits{std::move(tune_limits)}, tune_verbose{tune_verbose} {
 
     if (tune_timings <= 0) {
       throw std::domain_error("Parameter 'timings' must be > 0");

src/core/unit_tests/EspressoSystemStandAlone_test.cpp

@@ -75,6 +75,7 @@ namespace utf = boost::unit_test;
 #include <initializer_list>
 #include <limits>
 #include <memory>
+#include <optional>
 #include <stdexcept>
 #include <unordered_map>
 #include <utility>
@@ -228,6 +229,8 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory) {
 
     // set up P3M
     auto const prefactor = 2.;
+    auto const mesh_range = std::pair<std::optional<int>, std::optional<int>>{
+        std::nullopt, std::nullopt};
     auto p3m = P3MParameters{false,
                              0.0,
                              3.5,
@@ -238,7 +241,7 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory) {
                              1e-3};
     auto solver =
         new_p3m_handle<double, Arch::CPU, FFTBackendLegacy, FFTBuffersLegacy>(
-            std::move(p3m), prefactor, 1, false, true);
+            std::move(p3m), prefactor, 1, false, mesh_range, true);
     add_actor(comm, espresso::system, system.coulomb.impl->solver, solver,
               [&system]() { system.on_coulomb_change(); });
     BOOST_CHECK(not solver->is_gpu());
@@ -297,6 +300,8 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory) {
 
     // set up P3M
     auto const prefactor = 2.;
+    auto const mesh_range = std::pair<std::optional<int>, std::optional<int>>{
+        std::nullopt, std::nullopt};
     auto p3m = P3MParameters{false,
                              0.0,
                              3.5,
@@ -307,7 +312,7 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory) {
                              1e-3};
     auto solver =
         new_dp3m_handle<double, Arch::CPU, FFTBackendLegacy, FFTBuffersLegacy>(
-            std::move(p3m), prefactor, 1, false);
+            std::move(p3m), prefactor, 1, false, mesh_range);
     add_actor(comm, espresso::system, system.dipoles.impl->solver, solver,
               [&system]() { system.on_dipoles_change(); });
     BOOST_CHECK(not solver->is_gpu());

src/python/espressomd/electrostatics.py

@@ -209,6 +209,10 @@ class P3M(_P3MBase):
     tune : :obj:`bool`, optional
         Used to activate/deactivate the tuning method on activation.
         Defaults to ``True``.
+    tune_limits : (2,) array_like of :obj:`int`, optional
+        Lower and upper limits (inclusive) for the mesh size during tuning,
+        along the largest box dimension. Use ``None`` to not impose a limit.
+        Defaults to ``[None, None]``.
     timings : :obj:`int`
         Number of force calculations during tuning.
     verbose : :obj:`bool`, optional
@@ -266,6 +270,10 @@ class P3MGPU(_P3MBase):
         Defaults to ``True``.
     timings : :obj:`int`
         Number of force calculations during tuning.
+    tune_limits : (2,) array_like of :obj:`int`, optional
+        Lower and upper limits (inclusive) for the mesh size during tuning,
+        along the largest box dimension. Use ``None`` to not impose a limit.
+        Defaults to ``[None, None]``.
     verbose : :obj:`bool`, optional
         If ``False``, disable log output during tuning.
     check_neutrality : :obj:`bool`, optional

src/python/espressomd/magnetostatics.py

@@ -104,6 +104,9 @@ class DipolarP3M(MagnetostaticInteraction):
     tune : :obj:`bool`, optional
         Activate/deactivate the tuning method on activation
         (default is ``True``, i.e., activated).
+    tune_limits : (2,) array_like of :obj:`int`, optional
+        Lower and upper limits (inclusive) for the mesh size during tuning.
+        Use ``None`` to not impose a limit. Defaults to ``[None, None]``.
     timings : :obj:`int`
         Number of force calculations during tuning.
     single_precision : :obj:`bool`

src/script_interface/electrostatics/CoulombP3M.hpp

@@ -33,6 +33,7 @@
 #include "script_interface/get_value.hpp"
 
 #include <memory>
+#include <optional>
 #include <stdexcept>
 #include <string>
 #include <utility>
@@ -51,6 +52,7 @@ class CoulombP3M : public Actor<CoulombP3M<Architecture>, ::CoulombP3M> {
   bool m_tune_verbose;
   bool m_check_complex_residuals;
   bool m_single_precision;
+  std::pair<std::optional<int>, std::optional<int>> m_tune_limits;
 
 public:
   using Base = Actor<CoulombP3M<Architecture>, ::CoulombP3M>;
@@ -93,6 +95,15 @@ class CoulombP3M : public Actor<CoulombP3M<Architecture>, ::CoulombP3M> {
          [this]() { return m_tune_verbose; }},
         {"timings", AutoParameter::read_only,
          [this]() { return m_tune_timings; }},
+        {"tune_limits", AutoParameter::read_only,
+         [this]() {
+           auto const &[range_min, range_max] = m_tune_limits;
+           std::vector<Variant> retval = {
+               range_min ? Variant{*range_min} : Variant{None{}},
+               range_max ? Variant{*range_max} : Variant{None{}},
+           };
+           return retval;
+         }},
         {"tune", AutoParameter::read_only, [this]() { return m_tune; }},
         {"check_complex_residuals", AutoParameter::read_only,
          [this]() { return m_check_complex_residuals; }},
@@ -103,6 +114,33 @@ class CoulombP3M : public Actor<CoulombP3M<Architecture>, ::CoulombP3M> {
     m_tune = get_value<bool>(params, "tune");
     m_tune_timings = get_value<int>(params, "timings");
     m_tune_verbose = get_value<bool>(params, "verbose");
+    m_tune_limits = {std::nullopt, std::nullopt};
+    if (params.contains("tune_limits")) {
+      std::vector<Variant> range;
+      try {
+        auto const val = get_value<std::vector<int>>(params, "tune_limits");
+        assert(val.size() == 2u);
+        range.emplace_back(val[0u]);
+        range.emplace_back(val[1u]);
+      } catch (...) {
+        range = get_value<std::vector<Variant>>(params, "tune_limits");
+        assert(range.size() == 2u);
+      }
+      if (not is_none(range[0u])) {
+        m_tune_limits.first = get_value<int>(range[0u]);
+      }
+      if (not is_none(range[1u])) {
+        m_tune_limits.second = get_value<int>(range[1u]);
+      }
+      context()->parallel_try_catch([&]() {
+        if (m_tune_limits.first and *m_tune_limits.first <= 0) {
+          throw std::domain_error("P3M mesh tuning limits: mesh must be > 0");
+        }
+        if (m_tune_limits.second and *m_tune_limits.second <= 0) {
+          throw std::domain_error("P3M mesh tuning limits: mesh must be > 0");
+        }
+      });
+    }
     m_check_complex_residuals =
         get_value<bool>(params, "check_complex_residuals");
     auto const single_precision = get_value<bool>(params, "single_precision");
@@ -121,7 +159,7 @@ class CoulombP3M : public Actor<CoulombP3M<Architecture>, ::CoulombP3M> {
                                get_value<double>(params, "accuracy")};
       make_handle(single_precision, std::move(p3m),
                   get_value<double>(params, "prefactor"), m_tune_timings,
-                  m_tune_verbose, m_check_complex_residuals);
+                  m_tune_verbose, m_tune_limits, m_check_complex_residuals);
     });
     set_charge_neutrality_tolerance(params);
   }

src/script_interface/magnetostatics/DipolarP3M.hpp

@@ -33,6 +33,7 @@
 #include "script_interface/get_value.hpp"
 
 #include <memory>
+#include <optional>
 #include <stdexcept>
 #include <string>
 #include <utility>
@@ -47,6 +48,7 @@ namespace Dipoles {
 template <Arch Architecture>
 class DipolarP3M : public Actor<DipolarP3M<Architecture>, ::DipolarP3M> {
   int m_tune_timings;
+  std::pair<std::optional<int>, std::optional<int>> m_tune_limits;
   bool m_tune;
   bool m_tune_verbose;
 
@@ -90,6 +92,15 @@ class DipolarP3M : public Actor<DipolarP3M<Architecture>, ::DipolarP3M> {
          [this]() { return m_tune_verbose; }},
         {"timings", AutoParameter::read_only,
          [this]() { return m_tune_timings; }},
+        {"tune_limits", AutoParameter::read_only,
+         [this]() {
+           auto const &[range_min, range_max] = m_tune_limits;
+           std::vector<Variant> retval = {
+               range_min ? Variant{*range_min} : Variant{None{}},
+               range_max ? Variant{*range_max} : Variant{None{}},
+           };
+           return retval;
+         }},
         {"tune", AutoParameter::read_only, [this]() { return m_tune; }},
     });
   }
@@ -98,6 +109,33 @@ class DipolarP3M : public Actor<DipolarP3M<Architecture>, ::DipolarP3M> {
     m_tune = get_value<bool>(params, "tune");
     m_tune_timings = get_value<int>(params, "timings");
     m_tune_verbose = get_value<bool>(params, "verbose");
+    m_tune_limits = {std::nullopt, std::nullopt};
+    if (params.contains("tune_limits")) {
+      std::vector<Variant> range;
+      try {
+        auto const val = get_value<std::vector<int>>(params, "tune_limits");
+        assert(val.size() == 2u);
+        range.emplace_back(val[0u]);
+        range.emplace_back(val[1u]);
+      } catch (...) {
+        range = get_value<std::vector<Variant>>(params, "tune_limits");
+        assert(range.size() == 2u);
+      }
+      if (not is_none(range[0u])) {
+        m_tune_limits.first = get_value<int>(range[0u]);
+      }
+      if (not is_none(range[1u])) {
+        m_tune_limits.second = get_value<int>(range[1u]);
+      }
+      context()->parallel_try_catch([&]() {
+        if (m_tune_limits.first and *m_tune_limits.first <= 0) {
+          throw std::domain_error("P3M mesh tuning limits: mesh must be > 0");
+        }
+        if (m_tune_limits.second and *m_tune_limits.second <= 0) {
+          throw std::domain_error("P3M mesh tuning limits: mesh must be > 0");
+        }
+      });
+    }
     auto const single_precision = get_value<bool>(params, "single_precision");
     static_assert(Architecture == Arch::CPU, "GPU not implemented");
     context()->parallel_try_catch([&]() {
@@ -111,7 +149,7 @@ class DipolarP3M : public Actor<DipolarP3M<Architecture>, ::DipolarP3M> {
                                get_value<double>(params, "accuracy")};
       make_handle(single_precision, std::move(p3m),
                   get_value<double>(params, "prefactor"), m_tune_timings,
-                  m_tune_verbose);
+                  m_tune_verbose, m_tune_limits);
     });
   }