Simplify P3M C++ API (#4969)

Fixes #4962
Partial fix for #4963

Description of changes:
- disentangle `FFTBackend` methods: one per FFT transform operation, one per halo communication operation
- make the floating-point type and CPU/GPU architecture template parameters of P3M classes
- add support for P3M CPU and dipolar P3M CPU in single-precision mode

cmake/FindFFTW3.cmake

@@ -31,25 +31,27 @@ if(FFTW3_INCLUDE_DIR)
 endif(FFTW3_INCLUDE_DIR)
 
 find_path(FFTW3_INCLUDE_DIR fftw3.h)
-find_library(FFTW3_LIBRARIES NAMES fftw3)
+find_library(FFTW3_LIBRARIES fftw3)
+find_library(FFTW3F_LIBRARIES fftw3f)
 find_path(FFTW3_MPI_INCLUDE_DIR fftw3-mpi.h)
-find_library(FFTW3_MPI_LIBRARIES NAMES fftw3_mpi)
+find_library(FFTW3_MPI_LIBRARIES fftw3_mpi)
+find_library(FFTW3F_MPI_LIBRARIES fftw3f_mpi)
 
 # handle the QUIETLY and REQUIRED arguments and set FFTW3_FOUND to TRUE if all
 # listed variables are TRUE
 include(FindPackageHandleStandardArgs)
-find_package_handle_standard_args(FFTW3 DEFAULT_MSG FFTW3_LIBRARIES
+find_package_handle_standard_args(FFTW3 DEFAULT_MSG FFTW3_LIBRARIES FFTW3F_LIBRARIES
                                   FFTW3_INCLUDE_DIR)
 set(FPHSA_NAME_MISMATCHED 1)
-find_package_handle_standard_args(FFTW3_MPI DEFAULT_MSG FFTW3_MPI_LIBRARIES
+find_package_handle_standard_args(FFTW3_MPI DEFAULT_MSG FFTW3_MPI_LIBRARIES FFTW3F_MPI_LIBRARIES
                                   FFTW3_MPI_INCLUDE_DIR)
 unset(FPHSA_NAME_MISMATCHED)
 
-mark_as_advanced(FFTW3_LIBRARIES FFTW3_INCLUDE_DIR FFTW3_MPI_LIBRARIES FFTW3_MPI_INCLUDE_DIR)
+mark_as_advanced(FFTW3_LIBRARIES FFTW3F_LIBRARIES FFTW3_INCLUDE_DIR FFTW3_MPI_LIBRARIES FFTW3F_MPI_LIBRARIES FFTW3_MPI_INCLUDE_DIR)
 
 
 if(FFTW3_FOUND AND NOT TARGET FFTW3::FFTW3)
   add_library(FFTW3::FFTW3 INTERFACE IMPORTED)
   target_include_directories(FFTW3::FFTW3 INTERFACE "${FFTW3_INCLUDE_DIR}")
-  target_link_libraries(FFTW3::FFTW3 INTERFACE "${FFTW3_LIBRARIES}")
+  target_link_libraries(FFTW3::FFTW3 INTERFACE "${FFTW3_LIBRARIES}" "${FFTW3F_LIBRARIES}")
 endif()

src/core/electrostatics/CMakeLists.txt

@@ -24,6 +24,5 @@ target_sources(
           ${CMAKE_CURRENT_SOURCE_DIR}/mmm1d.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/mmm-modpsi.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/p3m.cpp
-          ${CMAKE_CURRENT_SOURCE_DIR}/p3m_gpu.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/scafacos_impl.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/specfunc.cpp)

src/core/electrostatics/coulomb.cpp

@@ -207,11 +207,6 @@ struct LongRangeForce {
   void operator()(std::shared_ptr<CoulombP3M> const &actor) const {
     actor->add_long_range_forces(m_particles);
   }
-#ifdef CUDA
-  void operator()(std::shared_ptr<CoulombP3MGPU> const &actor) const {
-    actor->add_long_range_forces(m_particles);
-  }
-#endif // CUDA
   void
   operator()(std::shared_ptr<ElectrostaticLayerCorrection> const &actor) const {
     actor->add_long_range_forces(m_particles);

src/core/electrostatics/coulomb.hpp

@@ -32,7 +32,6 @@
 #include "electrostatics/icc.hpp"
 #include "electrostatics/mmm1d.hpp"
 #include "electrostatics/p3m.hpp"
-#include "electrostatics/p3m_gpu.hpp"
 #include "electrostatics/reaction_field.hpp"
 #include "electrostatics/scafacos.hpp"
 
@@ -49,9 +48,6 @@ using ElectrostaticsActor =
     std::variant<std::shared_ptr<DebyeHueckel>,
 #ifdef P3M
                  std::shared_ptr<CoulombP3M>,
-#ifdef CUDA
-                 std::shared_ptr<CoulombP3MGPU>,
-#endif // CUDA
                  std::shared_ptr<ElectrostaticLayerCorrection>,
 #endif // P3M
                  std::shared_ptr<CoulombMMM1D>,

src/core/electrostatics/elc.cpp

@@ -27,7 +27,6 @@
 
 #include "electrostatics/coulomb.hpp"
 #include "electrostatics/p3m.hpp"
-#include "electrostatics/p3m_gpu.hpp"
 
 #include "BoxGeometry.hpp"
 #include "Particle.hpp"
@@ -1107,39 +1106,34 @@ ElectrostaticLayerCorrection::ElectrostaticLayerCorrection(
   adapt_solver();
 }
 
-static void p3m_assign_image_charge(elc_data const &elc, CoulombP3M &p3m,
-                                    double q, Utils::Vector3d const &pos) {
-  if (pos[2] < elc.space_layer) {
-    auto const q_eff = elc.delta_mid_bot * q;
-    p3m.assign_charge(q_eff, {pos[0], pos[1], -pos[2]});
-  }
-  if (pos[2] > (elc.box_h - elc.space_layer)) {
-    auto const q_eff = elc.delta_mid_top * q;
-    p3m.assign_charge(q_eff, {pos[0], pos[1], 2. * elc.box_h - pos[2]});
-  }
-}
-
 template <ChargeProtocol protocol, typename combined_ranges>
 void charge_assign(elc_data const &elc, CoulombP3M &solver,
                    combined_ranges const &p_q_pos_range) {
-  if (protocol == ChargeProtocol::BOTH or protocol == ChargeProtocol::IMAGE) {
-    solver.p3m.inter_weights.reset(solver.p3m.params.cao);
-  }
-  /* prepare local FFT mesh */
-  for (int i = 0; i < solver.p3m.local_mesh.size; i++)
-    solver.p3m.mesh.rs_scalar[i] = 0.;
+
+  solver.prepare_fft_mesh(protocol == ChargeProtocol::BOTH or
+                          protocol == ChargeProtocol::IMAGE);
 
   for (auto zipped : p_q_pos_range) {
     auto const p_q = boost::get<0>(zipped);
     auto const &p_pos = boost::get<1>(zipped);
     if (p_q != 0.) {
+      // assign real charges
       if (protocol == ChargeProtocol::BOTH or
           protocol == ChargeProtocol::REAL) {
-        solver.assign_charge(p_q, p_pos, solver.p3m.inter_weights);
+        solver.assign_charge(p_q, p_pos, false);
       }
+      // assign image charges
       if (protocol == ChargeProtocol::BOTH or
           protocol == ChargeProtocol::IMAGE) {
-        p3m_assign_image_charge(elc, solver, p_q, p_pos);
+        if (p_pos[2] < elc.space_layer) {
+          auto const q_eff = elc.delta_mid_bot * p_q;
+          solver.assign_charge(q_eff, {p_pos[0], p_pos[1], -p_pos[2]}, true);
+        }
+        if (p_pos[2] > (elc.box_h - elc.space_layer)) {
+          auto const q_eff = elc.delta_mid_top * p_q;
+          solver.assign_charge(
+              q_eff, {p_pos[0], p_pos[1], 2. * elc.box_h - p_pos[2]}, true);
+        }
       }
     }
   }
@@ -1149,7 +1143,9 @@ template <ChargeProtocol protocol, typename combined_range>
 void modify_p3m_sums(elc_data const &elc, CoulombP3M &solver,
                      combined_range const &p_q_pos_range) {
 
-  Utils::Vector3d node_sums{};
+  auto local_n = 0;
+  auto local_q2 = 0.0;
+  auto local_q = 0.0;
   for (auto zipped : p_q_pos_range) {
     auto const p_q = boost::get<0>(zipped);
     auto const &p_pos = boost::get<1>(zipped);
@@ -1158,33 +1154,35 @@ void modify_p3m_sums(elc_data const &elc, CoulombP3M &solver,
 
       if (protocol == ChargeProtocol::BOTH or
           protocol == ChargeProtocol::REAL) {
-        node_sums[0] += 1.;
-        node_sums[1] += Utils::sqr(p_q);
-        node_sums[2] += p_q;
+        local_n++;
+        local_q2 += Utils::sqr(p_q);
+        local_q += p_q;
       }
 
       if (protocol == ChargeProtocol::BOTH or
           protocol == ChargeProtocol::IMAGE) {
         if (p_z < elc.space_layer) {
-          node_sums[0] += 1.;
-          node_sums[1] += Utils::sqr(elc.delta_mid_bot * p_q);
-          node_sums[2] += elc.delta_mid_bot * p_q;
+          local_n++;
+          local_q2 += Utils::sqr(elc.delta_mid_bot * p_q);
+          local_q += elc.delta_mid_bot * p_q;
         }
 
         if (p_z > (elc.box_h - elc.space_layer)) {
-          node_sums[0] += 1.;
-          node_sums[1] += Utils::sqr(elc.delta_mid_top * p_q);
-          node_sums[2] += elc.delta_mid_top * p_q;
+          local_n++;
+          local_q2 += Utils::sqr(elc.delta_mid_top * p_q);
+          local_q += elc.delta_mid_top * p_q;
         }
       }
     }
   }
 
-  auto const tot_sums =
-      boost::mpi::all_reduce(comm_cart, node_sums, std::plus<>());
-  solver.p3m.sum_qpart = static_cast<int>(tot_sums[0] + 0.1);
-  solver.p3m.sum_q2 = tot_sums[1];
-  solver.p3m.square_sum_q = Utils::sqr(tot_sums[2]);
+  auto global_n = 0;
+  auto global_q2 = 0.;
+  auto global_q = 0.;
+  boost::mpi::all_reduce(comm_cart, local_n, global_n, std::plus<>());
+  boost::mpi::all_reduce(comm_cart, local_q2, global_q2, std::plus<>());
+  boost::mpi::all_reduce(comm_cart, local_q, global_q, std::plus<>());
+  solver.count_charged_particles_elc(global_n, global_q2, Utils::sqr(global_q));
 }
 
 double ElectrostaticLayerCorrection::long_range_energy(

src/core/electrostatics/elc.hpp

@@ -37,13 +37,13 @@
 #include "actor/traits.hpp"
 
 #include "electrostatics/p3m.hpp"
-#include "electrostatics/p3m_gpu.hpp"
 
 #include "BoxGeometry.hpp"
 #include "Particle.hpp"
 #include "ParticleRange.hpp"
 
 #include <utils/Vector.hpp>
+#include <utils/math/sqr.hpp>
 
 #include <memory>
 #include <optional>
@@ -166,11 +166,7 @@ struct elc_data {
 
 struct ElectrostaticLayerCorrection
     : public Coulomb::Actor<ElectrostaticLayerCorrection> {
-  using BaseSolver = std::variant<
-#ifdef CUDA
-      std::shared_ptr<CoulombP3MGPU>,
-#endif // CUDA
-      std::shared_ptr<CoulombP3M>>;
+  using BaseSolver = std::variant<std::shared_ptr<CoulombP3M>>;
 
   elc_data elc;
   BoxGeometry *m_box_geo;

src/core/electrostatics/icc.cpp

@@ -40,6 +40,7 @@
 #include "communication.hpp"
 #include "electrostatics/coulomb.hpp"
 #include "electrostatics/coulomb_inline.hpp"
+#include "electrostatics/p3m.hpp"
 #include "errorhandling.hpp"
 #include "integrators/Propagation.hpp"
 #include "system/System.hpp"
@@ -248,8 +249,10 @@ struct SanityChecksICC {
   template <typename T> void operator()(std::shared_ptr<T> const &) const {}
 #ifdef P3M
 #ifdef CUDA
-  [[noreturn]] void operator()(std::shared_ptr<CoulombP3MGPU> const &) const {
-    throw std::runtime_error("ICC does not work with P3MGPU");
+  void operator()(std::shared_ptr<CoulombP3M> const &p) const {
+    if (p->is_gpu()) {
+      throw std::runtime_error("ICC does not work with P3MGPU");
+    }
   }
 #endif // CUDA
   void