Integrate Exscalicos

CMakeLists.txt

@@ -104,6 +104,7 @@ option(ESPRESSO_BUILD_WITH_CUDA "Build with GPU support" OFF)
 option(ESPRESSO_BUILD_WITH_HDF5 "Build with HDF5 support" OFF)
 option(ESPRESSO_BUILD_TESTS "Enable tests" ON)
 option(ESPRESSO_BUILD_WITH_SCAFACOS "Build with ScaFaCoS support" OFF)
+option(ESPRESSO_BUILD_WITH_EXSCALICOS "Build with ExScaLiCoS support" OFF)
 option(ESPRESSO_BUILD_WITH_STOKESIAN_DYNAMICS "Build with Stokesian Dynamics"
        OFF)
 option(ESPRESSO_BUILD_WITH_NLOPT "Build with NLopt support" OFF)
@@ -742,6 +743,7 @@ if(NOT DEFINED Cabana_FOUND OR NOT ${Cabana_FOUND})
     GIT_REPOSITORY https://github.com/ECP-copa/Cabana.git
     GIT_TAG        e76c1a1 # 0.7.0 with patches
     PATCH_COMMAND  patch -p0 < ${CMAKE_CURRENT_SOURCE_DIR}/cmake/cabana.patch
+    OVERRIDE_FIND_PACKAGE
   )
   # cmake-format: on
   set(Cabana_REQUIRE_HEFFTE ${ESPRESSO_BUILD_WITH_FFTW} CACHE BOOL "")
@@ -790,6 +792,32 @@ if(ESPRESSO_BUILD_WITH_SCAFACOS)
   pkg_check_modules(SCAFACOS scafacos REQUIRED)
 endif()
 
+if(ESPRESSO_BUILD_WITH_EXSCALICOS)
+  if(NOT DEFINED Exscalicos_FOUND OR NOT ${Exscalicos_FOUND})
+    if(NOT EXISTS ${FETCHCONTENT_BASE_DIR}/exscalicos-src)
+      find_package(Exscalicos QUIET)
+    endif()
+    if(NOT DEFINED Exscalicos_FOUND OR NOT ${Exscalicos_FOUND})
+      # cmake-format: off
+      FetchContent_Declare(
+        exscalicos
+        GIT_REPOSITORY https://gitlab.jsc.fz-juelich.de/SLMS/kokkos/exscalicos.git
+        GIT_TAG        bcd2c238 # P3M version
+        PATCH_COMMAND  patch -p1 < ${CMAKE_CURRENT_SOURCE_DIR}/cmake/exscalicos.patch
+      )
+      # cmake-format: on
+      set(EXSCALICOS_BUILD_TESTS off CACHE BOOL "")
+      set(EXSCALICOS_CUDA off CACHE BOOL "")
+      set(EXSCALICOS_HIP off CACHE BOOL "")
+      set(BUILD_SHARED_LIBS OFF)
+      FetchContent_MakeAvailable(exscalicos)
+      set(BUILD_SHARED_LIBS ${ESPRESSO_BUILD_SHARED_LIBS_DEFAULT})
+      set(CMAKE_SHARED_LIBRARY_PREFIX "${ESPRESSO_SHARED_LIBRARY_PREFIX}")
+      add_library(Exscalicos::Exscalicos ALIAS Exscalicos)
+    endif()
+  endif()
+endif()
+
 if(ESPRESSO_BUILD_WITH_GSL)
   find_package(GSL REQUIRED)
 endif()

src/scafacos/CMakeLists.txt

@@ -21,8 +21,11 @@ add_library(espresso_scafacos SHARED src/Scafacos.cpp src/Coulomb.cpp
                                      src/Dipoles.cpp)
 add_library(espresso::scafacos ALIAS espresso_scafacos)
 espresso_set_common_target_properties(espresso_scafacos)
-target_link_libraries(espresso_scafacos PUBLIC MPI::MPI_CXX
-                      PRIVATE ${SCAFACOS_LDFLAGS} espresso::compiler_flags)
+target_link_libraries(
+  espresso_scafacos PUBLIC MPI::MPI_CXX
+  PRIVATE ${SCAFACOS_LDFLAGS}
+          $<$<BOOL:${WALBERLA_BUILD_WITH_EXSCALICOS}>:Exscalicos::Exscalicos>
+          espresso::compiler_flags)
 
 target_include_directories(
   espresso_scafacos

src/scafacos/include/scafacos/Coulomb.hpp

@@ -60,6 +60,8 @@ struct Coulomb : public Scafacos {
     if (m_delegate_near_field) {
       fcs_float field;
       fcs_compute_near_field(m_handle, dist, &field);
+      // field = exObj.shortRangeField(dist);
+      field = 0.0;
       return field;
     }
 
@@ -71,6 +73,8 @@ struct Coulomb : public Scafacos {
     if (m_delegate_near_field) {
       fcs_float potential;
       fcs_compute_near_potential(m_handle, dist, &potential);
+      // potential = exObj.shortRangePotential(dist);
+      potential = 0.0;
       return potential;
     }
 

src/scafacos/include/scafacos/Scafacos.hpp

@@ -20,6 +20,7 @@
 #ifndef ESPRESSO_SRC_SCAFACOS_SCAFACOS_HPP
 #define ESPRESSO_SRC_SCAFACOS_SCAFACOS_HPP
 
+#include "Exscalicos.hpp"
 #include <fcs.h>
 
 #include <mpi.h>
@@ -50,6 +51,18 @@ struct Scafacos {
   /** Handle from the library */
   FCS m_handle;
 
+#ifdef EXSCALICOS_USING_CUDA
+  Exscalicos<ExscalicosSolver::P3M, double, Kokkos::Cuda, Kokkos::SharedSpace>
+      exObj;
+  ExscalicosData<double, Kokkos::Cuda, Kokkos::SharedSpace> exData;
+#else
+  Exscalicos<ExscalicosSolver::P3M, double, Kokkos::OpenMP, Kokkos::HostSpace>
+      exObj;
+  ExscalicosData<double, Kokkos::OpenMP, Kokkos::HostSpace> exData;
+#endif
+
+  int setupExscalicos;
+
 private:
   /** The method name */
   std::string m_method_name;

src/scafacos/src/Coulomb.cpp

@@ -23,6 +23,8 @@
 
 #include <fcs.h>
 
+#include <algorithm>
+#include <chrono>
 #include <stdexcept>
 #include <string>
 #include <utility>
@@ -81,17 +83,169 @@ void Coulomb::run(std::vector<double> &charges, std::vector<double> &positions,
   fields.resize(3ul * n_part);
   potentials.resize(n_part);
 
-  tune(charges, positions);
+  std::cout << "RUN: allocating copies, size = " << n_part << std::endl;
+  std::vector<double> positions2(3ul * n_part);
+  // positions2.reserve(positions.size());
+  std::vector<double> charges2(n_part);
+  // charges2.reserve(charges.size());
+  std::vector<double> fields2(3ul * n_part);
+  // fields2.reserve(fields.size());
+  std::vector<double> potentials2(n_part);
+  // potentials2.reserve(potentials.size());
+
+  std::cout << "RUN: copying data" << std::endl;
+  std::copy(positions.begin(), positions.end(), positions2.begin());
+  std::copy(charges.begin(), charges.end(), charges2.begin());
+  std::copy(fields.begin(), fields.end(), fields2.begin());
+  std::copy(potentials.begin(), potentials.end(), potentials2.begin());
+
+  std::cout << "RUN: MPI parameters" << std::endl;
+  int rank;
+  MPI_Comm_rank(exData.globalComm(), &rank);
+  int nranks;
+  MPI_Comm_size(exData.globalComm(), &nranks);
+  std::cout << "RUN: assigning sizes for partitions" << std::endl;
+  exObj.setParameter(ExscalicosP3MShortRangeProcess(), exData, 0, nranks - 1);
+  exObj.setParameter(ExscalicosP3MLongRangeProcess(), exData, 0, nranks - 1);
+
+  // std::cout << rank << " running ExScaLiCoS::tune now - using " << nranks
+  //           << " ranks : " << std::endl;
+  if (setupExscalicos != 0) {
+    // auto rMax = r_cut();
+    // exObj.setRMax(rMax);
+    exObj.tune(exData, positions2, charges2);
+    setupExscalicos = 0;
+  }
+
+  std::cout << "RUN: running ScaFaCoS" << std::endl;
+  std::cout << rank << " running fcs_run now: " << std::endl;
+  // std::cout << "RUN: checking charges / positions content:" << std::endl;
+  // for (int i = 0; i < std::min(10, (int)charges.size()); ++i)
+  //   std::cout << rank << " " << i << " " << positions[3 * i + 0] << " "
+  //             << positions[3 * i + 1] << " " << positions[3 * i + 2] << " "
+  //             << charges[i] << std::endl;
+
   auto const size = static_cast<int>(n_part);
+
   handle_error(fcs_run(m_handle, size, positions.data(), charges.data(),
                        fields.data(), potentials.data()));
+  MPI_Barrier(exData.globalComm());
+  const auto startScafacos{std::chrono::steady_clock::now()};
+  handle_error(fcs_run(m_handle, size, positions.data(), charges.data(),
+                       fields.data(), potentials.data()));
+  const auto endScafacos{std::chrono::steady_clock::now()};
+  const std::chrono::duration<double> timeScafacos{endScafacos - startScafacos};
+  MPI_Barrier(exData.globalComm());
+
+  // std::cout << "RUN: checking charges2 / positions2 content:" << std::endl;
+  // for (int i = 0; i < std::min(10, (int)charges2.size()); ++i)
+  //   std::cout << rank << " " << i << " " << positions2[3 * i + 0] << " "
+  //             << positions2[3 * i + 1] << " " << positions2[3 * i + 2] << "
+  //             "
+  //             << charges2[i] << std::endl;
+
+  std::cout << "RUN: running ExScaLiCoS" << std::endl;
+  std::cout << rank << " running ExScaLiCoS::run now: " << std::endl;
+  MPI_Barrier(exData.globalComm());
+  const auto startExscalicos{std::chrono::steady_clock::now()};
+  exObj.run(exData, positions2, charges2, potentials2, fields2);
+  const auto endExscalicos{std::chrono::steady_clock::now()};
+  const std::chrono::duration<double> timeExscalicos{endExscalicos -
+                                                     startExscalicos};
+  MPI_Barrier(exData.globalComm());
+
+  double tScafacos = timeScafacos.count();
+  double tExscalicos = timeExscalicos.count();
+
+  MPI_Allreduce(&tScafacos, MPI_IN_PLACE, 1, MPI_DOUBLE, MPI_MAX,
+                exData.globalComm());
+  MPI_Allreduce(&tExscalicos, MPI_IN_PLACE, 1, MPI_DOUBLE, MPI_MAX,
+                exData.globalComm());
+
+  std::vector<double> diff = {0.0, 0.0, 0.0, 0.0, 0.0};
+  for (auto i = 0; i < charges2.size(); ++i) {
+    diff[0] +=
+        (fields[3 * i] - fields2[3 * i]) * (fields[3 * i] - fields2[3 * i]);
+    diff[1] += (fields[3 * i + 1] - fields2[3 * i + 1]) *
+               (fields[3 * i + 1] - fields2[3 * i + 1]);
+    diff[2] += (fields[3 * i + 2] - fields2[3 * i + 2]) *
+               (fields[3 * i + 2] - fields2[3 * i + 2]);
+    diff[3] +=
+        (potentials[i] - potentials2[i]) * (potentials[i] - potentials2[i]);
+  }
+  diff[4] = (double)(charges.size());
+
+  MPI_Allreduce(diff.data(), MPI_IN_PLACE, 5, MPI_DOUBLE, MPI_SUM,
+                exData.globalComm());
+
+  diff[0] /= diff[4];
+  diff[1] /= diff[4];
+  diff[2] /= diff[4];
+  diff[3] /= diff[4];
+
+  diff[0] = std::sqrt(diff[0]);
+  diff[1] = std::sqrt(diff[1]);
+  diff[2] = std::sqrt(diff[2]);
+  diff[3] = std::sqrt(diff[3]);
+
+  if (rank == 0) {
+    std::cout << "Comparision ScaFaCoS - ExScaLiCoS:\n deltaP = " << diff[3]
+              << '\n'
+              << "deltaf = " << diff[0] << " " << diff[1] << " " << diff[2]
+              << '\n'
+              << "tScaFaCoS: " << tScafacos << "\n"
+              << "tExScaLiCoS: " << tExscalicos << "\n"
+              << std::endl;
+  }
+
+  /*
+  if (rank == 0) {
+    std::cout << "Comparison results (fields / potentials): " << std::endl;
+    for (int i = 0; i < potentials.size(); ++i) {
+      std::cout << " ( " << potentials[i] << " , " << potentials2[i] << " ) "
+                << " ( " << fields[3 * i + 0] << " , " << fields2[3 * i + 0]
+                << " ) "
+                << " ( " << fields[3 * i + 1] << " , " << fields2[3 * i + 1]
+                << " ) "
+                << " ( " << fields[3 * i + 2] << " , " << fields2[3 * i + 2]
+                << " ) " << std::endl;
+    }
+  }
+  */
+
+  fields = fields2;
+  potentials = potentials2;
 }
 
 void Coulomb::tune(std::vector<double> &charges,
                    std::vector<double> &positions) {
   auto const n_part = charges.size();
   auto const size = static_cast<int>(n_part);
+
   handle_error(fcs_tune(m_handle, size, positions.data(), charges.data()));
+
+  std::cout << "RUN: allocating copies, size = " << n_part << std::endl;
+  std::vector<double> positions2(3ul * n_part);
+  // positions2.reserve(positions.size());
+  std::vector<double> charges2(n_part);
+  // charges2.reserve(charges.size());
+
+  std::cout << "RUN: copying data" << std::endl;
+  std::copy(positions.begin(), positions.end(), positions2.begin());
+  std::copy(charges.begin(), charges.end(), charges2.begin());
+
+  std::cout << "RUN: MPI parameters" << std::endl;
+  int rank;
+  MPI_Comm_rank(exData.globalComm(), &rank);
+  int nranks;
+  MPI_Comm_size(exData.globalComm(), &nranks);
+  std::cout << "RUN: assigning sizes for partitions" << std::endl;
+  exObj.setParameter(ExscalicosP3MShortRangeProcess(), exData, 0, nranks - 1);
+  exObj.setParameter(ExscalicosP3MLongRangeProcess(), exData, 0, nranks - 1);
+
+  // exObj.tune(exData, charges2, positions2);
+
+  std::cout << "FINISHED tuning" << std::endl;
 }
 
 } // namespace Scafacos

src/scafacos/src/Scafacos.cpp

@@ -71,16 +71,24 @@ std::vector<std::string> Scafacos::available_methods() {
 }
 
 Scafacos::Scafacos(MPI_Comm comm, std::string method, std::string parameters)
-    : m_method_name{std::move(method)}, m_parameters{std::move(parameters)} {
+    : exObj(comm), m_method_name{std::move(method)},
+      m_parameters{std::move(parameters)} {
 
   handle_error(fcs_init(&m_handle, m_method_name.c_str(), comm));
 
   fcs_set_resort(m_handle, 0);
 
   handle_error(fcs_parser(m_handle, m_parameters.c_str(), 0));
+
+  exData.globalComm(comm);
+  exObj.create(exData);
+  setupExscalicos = 1;
 }
 
-Scafacos::~Scafacos() { fcs_destroy(m_handle); }
+Scafacos::~Scafacos() {
+  fcs_destroy(m_handle);
+  exObj.destroy(exData);
+}
 
 void Scafacos::set_runtime_parameters(double const *box_l,
                                       int const *periodicity,
@@ -93,6 +101,27 @@ void Scafacos::set_runtime_parameters(double const *box_l,
   double off[3] = {0., 0., 0.};
   handle_error(fcs_set_common(m_handle, near_field_flag, boxa, boxb, boxc, off,
                               periodicity, total_particles));
+
+  handle_error(fcs_set_common(m_handle, 1, boxa, boxb, boxc, off, periodicity,
+                              total_particles));
+
+  exData.sysSize(0, box_l[0]);
+  exData.sysSize(1, box_l[1]);
+  exData.sysSize(2, box_l[2]);
+
+  exData.periodicity(0, periodicity[0]);
+  exData.periodicity(1, periodicity[1]);
+  exData.periodicity(2, periodicity[2]);
+
+  exData.nSysParticles((long)total_particles);
+
+  int tol_type;
+  double tolerance;
+  fcs_get_tolerance(m_handle, &tol_type, &tolerance);
+  exData.accuracy(tolerance);
+
+  exObj.setNearFieldFlag(fcs_get_near_field_flag(m_handle) != 0);
+  exObj.setNearFieldFlag(true);
 }
 
 } // namespace Scafacos

testsuite/python/CMakeLists.txt

@@ -233,6 +233,11 @@ python_test(FILE cutoffs.py MAX_NUM_PROC 4)
 python_test(FILE cutoffs.py MAX_NUM_PROC 1 SUFFIX 1_core)
 python_test(FILE constraint_shape_based.py MAX_NUM_PROC 2)
 python_test(FILE coulomb_cloud_wall.py MAX_NUM_PROC 4 GPU_SLOTS 3)
+if(ESPRESSO_BUILD_WITH_EXSCALICOS)
+  python_test(FILE coulomb_cloud_wall_hp.py MAX_NUM_PROC 4 GPU_SLOTS 3)
+  python_test(FILE coulomb_dna_water.py MAX_NUM_PROC 4 GPU_SLOTS 3)
+  python_test(FILE coulomb_silica_melt.py MAX_NUM_PROC 4 GPU_SLOTS 3)
+endif()
 python_test(FILE coulomb_tuning.py MAX_NUM_PROC 4 GPU_SLOTS 3 LABELS long)
 python_test(FILE accumulator_correlator.py MAX_NUM_PROC 4)
 python_test(FILE accumulator_mean_variance.py MAX_NUM_PROC 4)

testsuite/python/coulomb_cloud_wall.py

@@ -39,7 +39,7 @@ class CoulombCloudWall(ut.TestCase):
         "coulomb_cloud_wall_system.data"))
 
     tolerance = 1E-3
-    p3m_params = {'r_cut': 1.001, 'accuracy': 1e-3,
+    p3m_params = {'r_cut': 1.001, 'accuracy': 1e-4,
                   'mesh': [64, 64, 64], 'cao': 7, 'alpha': 2.70746}
 
     # Reference energy from P3M
@@ -101,14 +101,28 @@ def test_scafacos_p2nfft(self):
         self.system.electrostatics.solver = espressomd.electrostatics.Scafacos(
             prefactor=2.8,
             method_name="p2nfft",
-            method_params={"p2nfft_r_cut": 1.001, "tolerance_field": 1E-4})
+            method_params={"p2nfft_r_cut": 1.001, "tolerance_field": 1E-5})
         self.system.integrator.run(0)
         self.compare(
             "scafacos_p2nfft",
             prefactor=2.8,
             force_tol=1e-3,
             energy_tol=1e-3)
 
+    @utx.skipIfMissingFeatures(["SCAFACOS"])
+    @utx.skipIfMissingScafacosMethod("p3m")
+    def test_scafacos_p3m(self):
+        self.system.electrostatics.solver = espressomd.electrostatics.Scafacos(
+            prefactor=2.8,
+            method_name="p3m",
+            method_params={"p3m_cao": 7, "tolerance_field": 1E-5})
+        self.system.integrator.run(0)
+        self.compare(
+            "scafacos_p3m",
+            prefactor=2.8,
+            force_tol=1e-3,
+            energy_tol=1e-3)
+
     def test_zz_deactivation(self):
         # Is the energy and force 0, if no methods active
         self.assertEqual(self.system.analysis.energy()["total"], 0.0)