Enable two ScaFaCoS instances (espressomd#4036)

Fixes espressomd#2636

Description of changes:
- create two global variables: one for ScaFaCoS electrostatics and one for ScaFaCoS magnetostatics
- ScaFaCoS can now be used for electrostatics+magnetostatics simultaneously (see test in 9b65a90)
- hide the ScaFaCoS context (i.e. flattened position/charge/dipole arrays) in a `ScafacosContext` struct
- make this struct available from outside only through a pointer to `ScafacosContextBase*` (type erasure)
- make the ScaFaCoS short-range force/energy kernels inline
- fix a bug in the `dipolar_mdlc_p3m_scafacos_p2nfft.py` testcase
- fix a bug in the Cython interface of the `Scafacos` actor that prevented returning some parameters from the core

doc/sphinx/electrostatics.rst

@@ -386,6 +386,4 @@ for an accuracy of :math:`10^{-3}`::
 
 For details of the various methods and their parameters please refer to
 the ScaFaCoS manual. To use this feature, ScaFaCoS has to be built as a
-shared library. ScaFaCoS can be used only once, either for Coulomb or for
-dipolar interactions.
-
+shared library.

doc/sphinx/installation.rst

@@ -303,7 +303,7 @@ General features
        :ref:`Magnetostatics / Dipolar interactions`
        :ref:`Electrostatics`
 
--  ``SCAFACOS_DIPOLES``
+-  ``SCAFACOS_DIPOLES`` This activates magnetostatics methods of ScaFaCoS.
 
 -  ``ROTATION`` Switch on rotational degrees of freedom for the particles, as well as
    the corresponding quaternion integrator.

doc/sphinx/magnetostatics.rst

@@ -211,6 +211,4 @@ root mean square accuracy for the magnetic field.
 
 For details of the various methods and their parameters please refer to
 the ScaFaCoS manual. To use this feature, ScaFaCoS has to be built as a
-shared library. ScaFaCoS can be used only once, either for Coulomb or for
-dipolar interactions.
-
+shared library.

src/core/electrostatics_magnetostatics/CMakeLists.txt

@@ -14,6 +14,7 @@ target_sources(
           ${CMAKE_CURRENT_SOURCE_DIR}/p3m-dipolar.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/p3m_gpu.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/scafacos.cpp
+          ${CMAKE_CURRENT_SOURCE_DIR}/ScafacosContext.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/fft.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/coulomb.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/dipole.cpp

src/core/electrostatics_magnetostatics/ScafacosContext.cpp

@@ -0,0 +1,247 @@
+/*
+ * Copyright (C) 2010-2020 The ESPResSo project
+ * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
+ *   Max-Planck-Institute for Polymer Research, Theory Group
+ *
+ * This file is part of ESPResSo.
+ *
+ * ESPResSo is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * ESPResSo is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+/** @file
+ *  Provide a C-like interface for ScaFaCoS.
+ */
+
+#include "config.hpp"
+
+#if defined(SCAFACOS)
+
+#include "electrostatics_magnetostatics/ScafacosContext.hpp"
+
+#include "cells.hpp"
+#include "communication.hpp"
+#include "electrostatics_magnetostatics/coulomb.hpp"
+#include "electrostatics_magnetostatics/dipole.hpp"
+#include "electrostatics_magnetostatics/scafacos.hpp"
+#include "grid.hpp"
+#include "integrate.hpp"
+#include "particle_data.hpp"
+#include "tuning.hpp"
+
+#include <utils/Vector.hpp>
+
+#include <boost/mpi/collectives.hpp>
+#include <boost/range/algorithm/min_element.hpp>
+#include <boost/range/numeric.hpp>
+
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+#include <functional>
+#include <limits>
+#include <list>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
+namespace Scafacos {
+
+std::string ScafacosContext::get_method_and_parameters() {
+  std::string representation = get_method() + " " + get_parameters();
+  std::replace(representation.begin(), representation.end(), ',', ' ');
+  return representation;
+}
+
+void ScafacosContext::update_system_params() {
+  int per[3] = {box_geo.periodic(0) != 0, box_geo.periodic(1) != 0,
+                box_geo.periodic(2) != 0};
+
+  auto const n_part = boost::mpi::all_reduce(
+      comm_cart, cell_structure.local_particles().size(), std::plus<>());
+  set_common_parameters(box_geo.length().data(), per, n_part);
+}
+
+void ScafacosContextCoulomb::update_particle_data() {
+  positions.clear();
+  charges.clear();
+
+  for (auto const &p : cell_structure.local_particles()) {
+    auto const pos = folded_position(p.r.p, box_geo);
+    positions.push_back(pos[0]);
+    positions.push_back(pos[1]);
+    positions.push_back(pos[2]);
+    charges.push_back(p.p.q);
+  }
+}
+
+#ifdef SCAFACOS_DIPOLES
+void ScafacosContextDipoles::update_particle_data() {
+  positions.clear();
+  dipoles.clear();
+
+  for (auto const &p : cell_structure.local_particles()) {
+    auto const pos = folded_position(p.r.p, box_geo);
+    positions.push_back(pos[0]);
+    positions.push_back(pos[1]);
+    positions.push_back(pos[2]);
+    auto const dip = p.calc_dip();
+    dipoles.push_back(dip[0]);
+    dipoles.push_back(dip[1]);
+    dipoles.push_back(dip[2]);
+  }
+}
+#endif
+
+void ScafacosContextCoulomb::update_particle_forces() const {
+  if (positions.empty())
+    return;
+
+  int it = 0;
+  for (auto &p : cell_structure.local_particles()) {
+    p.f.f += coulomb.prefactor * p.p.q *
+             Utils::Vector3d(Utils::Span<const double>(&(fields[it]), 3));
+    it += 3;
+  }
+
+  /* Check that the particle number did not change */
+  assert(it == fields.size());
+}
+
+#ifdef SCAFACOS_DIPOLES
+void ScafacosContextDipoles::update_particle_forces() const {
+  if (positions.empty())
+    return;
+
+  int it = 0;
+  for (auto &p : cell_structure.local_particles()) {
+    // Indices
+    // 3 "potential" values per particles (see below)
+    int const it_p = 3 * it;
+    // 6 "field" values per particles (see below)
+    int const it_f = 6 * it;
+
+    // The scafacos term "potential" here in fact refers to the magnetic
+    // field. So, the torques are given by m \times B
+    auto const dip = p.calc_dip();
+    auto const t = vector_product(
+        dip,
+        Utils::Vector3d(Utils::Span<const double>(&(potentials[it_p]), 3)));
+    // The force is given by G m, where G is a matrix
+    // which comes from the "fields" output of scafacos like this
+    // 0 1 2
+    // 1 3 4
+    // 2 4 5
+    // where the numbers refer to indices in the "field" output from scafacos
+    auto const G = Utils::Vector<Utils::Vector3d, 3>{
+        {fields[it_f + 0], fields[it_f + 1], fields[it_f + 2]},
+        {fields[it_f + 1], fields[it_f + 3], fields[it_f + 4]},
+        {fields[it_f + 2], fields[it_f + 4], fields[it_f + 5]}};
+    auto const f = G * dip;
+
+    // Add to particles
+    p.f.f += dipole.prefactor * f;
+    p.f.torque += dipole.prefactor * t;
+    it++;
+  }
+
+  /* Check that the particle number did not change */
+  assert(it == positions.size() / 3);
+}
+#endif
+
+double ScafacosContextCoulomb::long_range_energy() {
+  update_particle_data();
+  run(charges, positions, fields, potentials);
+  return 0.5 * coulomb.prefactor *
+         boost::inner_product(charges, potentials, 0.0);
+}
+
+#ifdef SCAFACOS_DIPOLES
+double ScafacosContextDipoles::long_range_energy() {
+  update_particle_data();
+  run_dipolar(dipoles, positions, fields, potentials);
+  return -0.5 * dipole.prefactor *
+         boost::inner_product(dipoles, potentials, 0.0);
+}
+#endif
+
+static void set_r_cut_and_tune_local_worker(double r_cut) {
+  set_r_cut_and_tune_local(r_cut);
+}
+
+REGISTER_CALLBACK(set_r_cut_and_tune_local_worker)
+
+/** Determine runtime for a specific cutoff */
+static double time_r_cut(double r_cut) {
+  /* Set cutoff to time */
+  mpi_call_all(set_r_cut_and_tune_local_worker, r_cut);
+
+  return time_force_calc(10);
+}
+
+/** Determine the optimal cutoff by bisection */
+static void tune_r_cut() {
+  const double tune_limit = 1e-3;
+
+  auto const min_box_l = *boost::min_element(box_geo.length());
+  auto const min_local_box_l = *boost::min_element(local_geo.length());
+
+  /* scafacos p3m and Ewald do not accept r_cut 0 for no good reason */
+  double r_min = 1.0;
+  double r_max = std::min(min_local_box_l, min_box_l / 2.0) - skin;
+  double t_min = 0;
+  double t_max = std::numeric_limits<double>::max();
+
+  /* Run bisection */
+  while (std::fabs(r_min - r_max) > tune_limit) {
+    const double dr = 0.5 * (r_max - r_min);
+    const double t_mid = time_r_cut(r_min + dr);
+    t_min = time_r_cut(r_min);
+    t_max = time_r_cut(r_max);
+
+    if (t_min <= 0.0 or t_max <= 0.0) {
+      break;
+    }
+
+    if (t_mid > t_min) {
+      r_max = r_min += dr;
+    } else {
+      r_min += dr;
+    }
+  }
+}
+
+void ScafacosContextCoulomb::tune() {
+  update_particle_data();
+
+  // Check whether we have to do a bisection for the short range cutoff
+  // Check if there is a user-supplied cutoff
+  if (has_near && r_cut() <= 0.0) {
+    // run tuning on the head node
+    if (this_node == 0) {
+      tune_r_cut();
+    }
+  } else {
+    // ESPResSo is not affected by a short range cutoff -> tune in parallel
+    Scafacos::tune(charges, positions);
+  }
+}
+
+void ScafacosContextCoulomb::set_r_cut_and_tune_local(double r_cut) {
+  update_particle_data();
+  set_r_cut(r_cut);
+  Scafacos::tune(charges, positions);
+}
+
+} // namespace Scafacos
+#endif /* SCAFACOS */

src/core/electrostatics_magnetostatics/ScafacosContext.hpp

@@ -0,0 +1,107 @@
+/*
+ * Copyright (C) 2010-2020 The ESPResSo project
+ * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
+ *   Max-Planck-Institute for Polymer Research, Theory Group
+ *
+ * This file is part of ESPResSo.
+ *
+ * ESPResSo is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * ESPResSo is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+#ifndef SRC_CORE_ELECTROSTATICS_MAGNETOSTATICS_SCAFACOSCONTEXT_HPP
+#define SRC_CORE_ELECTROSTATICS_MAGNETOSTATICS_SCAFACOSCONTEXT_HPP
+/**
+ * @file
+ * @ref Scafacos::ScafacosContextBase implements the interface of the
+ * ScaFaCoS bridge. It is further derived for the coulombic and dipolar
+ * versions of ScaFaCoS.
+ */
+
+#include "config.hpp"
+
+#if defined(SCAFACOS)
+
+#include "Scafacos.hpp"
+#include "electrostatics_magnetostatics/ScafacosContextBase.hpp"
+
+#include <utils/Vector.hpp>
+
+#include <string>
+#include <vector>
+
+#if defined(SCAFACOS_DIPOLES) && !defined(FCS_ENABLE_DIPOLES)
+#error                                                                         \
+    "SCAFACOS_DIPOLES requires dipoles support in scafacos library (FCS_ENABLE_DIPOLES)."
+#endif
+
+namespace Scafacos {
+
+/** Encapsulation for the particle data needed by ScaFaCoS */
+struct ScafacosContext : ScafacosContextBase, Scafacos {
+  using Scafacos::Scafacos;
+  using ScafacosContextBase::ScafacosContextBase;
+  ~ScafacosContext() override = default;
+  void update_system_params() override;
+  double get_pair_force(double dist) const override {
+    return Scafacos::pair_force(dist);
+  }
+  double get_pair_energy(double dist) const override {
+    return Scafacos::pair_energy(dist);
+  }
+  std::string get_method_and_parameters() override;
+  double get_r_cut() const override { return Scafacos::r_cut(); }
+
+protected:
+  /** Outputs */
+  std::vector<double> fields, potentials;
+};
+
+struct ScafacosContextCoulomb : ScafacosContext {
+  using ScafacosContext::ScafacosContext;
+  void update_particle_data() override;
+  void update_particle_forces() const override;
+  double long_range_energy() override;
+  void add_long_range_force() override {
+    update_particle_data();
+    run(charges, positions, fields, potentials);
+    update_particle_forces();
+  }
+  void tune();
+  void set_r_cut_and_tune_local(double r_cut);
+
+private:
+  /** Inputs */
+  std::vector<double> positions, charges;
+};
+
+#ifdef SCAFACOS_DIPOLES
+struct ScafacosContextDipoles : ScafacosContext {
+  using ScafacosContext::ScafacosContext;
+  void update_particle_data() override;
+  void update_particle_forces() const override;
+  double long_range_energy() override;
+  void add_long_range_force() override {
+    update_particle_data();
+    run_dipolar(dipoles, positions, fields, potentials);
+    update_particle_forces();
+  }
+
+private:
+  /** Inputs */
+  std::vector<double> positions, dipoles;
+};
+#endif
+
+} // namespace Scafacos
+#endif // SCAFACOS
+#endif // SRC_CORE_ELECTROSTATICS_MAGNETOSTATICS_SCAFACOSCONTEXT_HPP