IBM volume conservation maintenance (#4286)

Fixes #4284

Description of changes:
- use a `std::vector` for memory management
- remove `IBM_MAX_NUM` macro, now that any `softID` is valid
- don't crash the simulation when measuring the system energy with IBM bonds
- write a testcase for the IBM volume conservation (partial fix for #2071)

src/config/config.hpp

@@ -100,9 +100,4 @@
 #define MAX_OBJECTS_IN_FLUID 10000
 #endif
 
-/** Maximal number of soft particles per immersed boundary */
-#ifndef IBM_MAX_NUM
-#define IBM_MAX_NUM 1000
-#endif
-
 #endif

src/core/energy_inline.hpp

@@ -257,6 +257,12 @@ calc_bonded_energy(Bonded_IA_Parameters const &iaparams, Particle const &p1,
     if (auto const *iap = boost::get<TabulatedAngleBond>(&iaparams)) {
       return iap->energy(p1.r.p, p2->r.p, p3->r.p);
     }
+    if (auto const *iap = boost::get<IBMTriel>(&iaparams)) {
+      runtimeWarningMsg() << "Unsupported bond type " +
+                                 std::to_string(iaparams.which()) +
+                                 " in energy calculation.";
+      return 0.;
+    }
     throw BondUnknownTypeError();
   } // 2 partners
   if (n_partners == 3) {
@@ -266,6 +272,12 @@ calc_bonded_energy(Bonded_IA_Parameters const &iaparams, Particle const &p1,
     if (auto const *iap = boost::get<TabulatedDihedralBond>(&iaparams)) {
       return iap->energy(p2->r.p, p1.r.p, p3->r.p, p4->r.p);
     }
+    if (auto const *iap = boost::get<IBMTribend>(&iaparams)) {
+      runtimeWarningMsg() << "Unsupported bond type " +
+                                 std::to_string(iaparams.which()) +
+                                 " in energy calculation.";
+      return 0.;
+    }
     throw BondUnknownTypeError();
   } // 3 partners
   if (n_partners == 0) {

src/core/immersed_boundary/ImmersedBoundaries.cpp

@@ -32,10 +32,11 @@
 #include <utils/Vector.hpp>
 #include <utils/constants.hpp>
 
-#include <mpi.h>
+#include <boost/mpi/collectives/all_reduce.hpp>
 
-#include <cstdio>
+#include <functional>
 #include <utility>
+#include <vector>
 
 /** Calculate volumes, volume force and add it to each virtual particle. */
 void ImmersedBoundaries::volume_conservation(CellStructure &cs) {
@@ -69,7 +70,7 @@ void ImmersedBoundaries::init_volume_conservation(CellStructure &cs) {
         // accidentally during the integration. Then we must not reset the
         // reference
         BoundariesFound = true;
-        if (v->volRef == 0) {
+        if (v->volRef == 0.) {
           v->volRef = VolumesCurrent[v->softID];
         }
       }
@@ -98,7 +99,7 @@ void ImmersedBoundaries::calc_volumes(CellStructure &cs) {
     return;
 
   // Partial volumes for each soft particle, to be summed up
-  std::vector<double> tempVol(IBM_MAX_NUM);
+  std::vector<double> tempVol(VolumesCurrent.size());
 
   // Loop over all particles on local node
   cs.bond_loop(
@@ -143,12 +144,10 @@ void ImmersedBoundaries::calc_volumes(CellStructure &cs) {
         return false;
       });
 
-  for (int i = 0; i < IBM_MAX_NUM; i++)
-    VolumesCurrent[i] = 0;
-
   // Sum up and communicate
-  MPI_Allreduce(tempVol.data(), VolumesCurrent.data(), IBM_MAX_NUM, MPI_DOUBLE,
-                MPI_SUM, comm_cart);
+  boost::mpi::all_reduce(comm_cart, tempVol.data(),
+                         static_cast<int>(tempVol.size()),
+                         VolumesCurrent.data(), std::plus<double>());
 }
 
 /** Calculate and add the volume force to each node */

src/core/immersed_boundary/ImmersedBoundaries.hpp

@@ -23,15 +23,29 @@
 
 #include "CellStructure.hpp"
 
+#include <cassert>
+#include <cstddef>
 #include <vector>
 
 class ImmersedBoundaries {
 public:
   ImmersedBoundaries() : VolumeInitDone(false), BoundariesFound(false) {
-    VolumesCurrent.resize(IBM_MAX_NUM);
+    VolumesCurrent.resize(10);
   }
   void init_volume_conservation(CellStructure &cs);
   void volume_conservation(CellStructure &cs);
+  void register_softID(int softID) {
+    assert(softID >= 0);
+    auto const new_size = static_cast<std::size_t>(softID) + 1;
+    if (new_size > VolumesCurrent.size()) {
+      VolumesCurrent.resize(new_size);
+    }
+  }
+  double get_current_volume(int softID) const {
+    assert(softID >= 0);
+    assert(softID < VolumesCurrent.size());
+    return VolumesCurrent[softID];
+  }
 
 private:
   void calc_volumes(CellStructure &cs);

src/core/immersed_boundary/ibm_volcons.cpp

@@ -18,25 +18,31 @@
  */
 
 #include "ibm_volcons.hpp"
-#include "errorhandling.hpp"
+
+#include "communication.hpp"
+#include "immersed_boundaries.hpp"
+
+#include <cassert>
+#include <stdexcept>
+
+void mpi_set_n_ibm_volcons_bonds_local(int softID) {
+  immersed_boundaries.register_softID(softID);
+}
+
+REGISTER_CALLBACK(mpi_set_n_ibm_volcons_bonds_local)
+
+void mpi_set_n_ibm_volcons_bonds(int softID) {
+  mpi_call_all(mpi_set_n_ibm_volcons_bonds_local, softID);
+}
 
 /** Set parameters of volume conservation */
 IBMVolCons::IBMVolCons(const int softID, const double kappaV) {
-  // Specific stuff
-  if (softID > IBM_MAX_NUM) {
-    runtimeErrorMsg() << "Error: softID (" << softID
-                      << ") is larger than IBM_MAX_NUM (" << IBM_MAX_NUM << ")";
-  }
-  if (softID < 0) {
-    runtimeErrorMsg() << "Error: softID (" << softID
-                      << ") must be non-negative";
-  }
-
   this->softID = softID;
   this->kappaV = kappaV;
-  volRef = 0;
   // NOTE: We cannot compute the reference volume here because not all
   // interactions are setup and thus we do not know which triangles belong to
-  // this softID. Calculate it later in the init function of \ref
-  // ImmersedBoundaries
+  // this softID. Calculate it later in the init function of
+  // \ref ImmersedBoundaries::init_volume_conservation()
+  volRef = 0.;
+  mpi_set_n_ibm_volcons_bonds(softID);
 }

src/python/espressomd/interactions.pxd

@@ -528,7 +528,7 @@ cdef extern from "thermostat.hpp":
 
 cdef extern from "immersed_boundary/ImmersedBoundaries.hpp":
     cppclass ImmersedBoundaries:
-        pass
+        double get_current_volume(int softID)
 
 cdef extern from "immersed_boundaries.hpp":
     extern ImmersedBoundaries immersed_boundaries

src/python/espressomd/interactions.pyx

@@ -21,7 +21,7 @@ from cython.operator cimport dereference
 import collections
 
 include "myconfig.pxi"
-from .utils import requires_experimental_features, is_valid_type
+from .utils import is_valid_type
 from .utils cimport check_type_or_throw_except
 
 
@@ -2942,7 +2942,6 @@ class IBM_Tribend(BondedInteraction):
             self._params["ind4"], self._params["kb"], flat))
 
 
-@requires_experimental_features("No test coverage")
 class IBM_VolCons(BondedInteraction):
 
     """
@@ -2954,7 +2953,8 @@ class IBM_VolCons(BondedInteraction):
     ----------
     softID : :obj:`int`
         Used to identify the object to which this bond belongs. Each object
-        (cell) needs its own ID
+        (cell) needs its own ID. For performance reasons, it is best to
+        start from ``softID=0`` and increment by 1 for each subsequent bond.
     kappaV : :obj:`float`
         Modulus for volume force
 
@@ -2988,6 +2988,14 @@ class IBM_VolCons(BondedInteraction):
         set_bonded_ia_params(self._bond_id, CoreIBMVolCons(
             self._params["softID"], self._params["kappaV"]))
 
+    def current_volume(self):
+        """
+        Query the current volume of the soft object associated to this bond.
+        The volume is initialized once all :class:`IBM_Triel` bonds have
+        been added and the forces have been recalculated.
+        """
+        return immersed_boundaries.get_current_volume(self._params['softID'])
+
 
 class OifGlobalForces(BondedInteraction):
 

testsuite/python/CMakeLists.txt

@@ -142,6 +142,7 @@ python_test(FILE domain_decomposition.py MAX_NUM_PROC 4)
 python_test(FILE integrator_npt.py MAX_NUM_PROC 4)
 python_test(FILE integrator_npt_stats.py MAX_NUM_PROC 4 LABELS long)
 python_test(FILE integrator_steepest_descent.py MAX_NUM_PROC 4)
+python_test(FILE ibm.py MAX_NUM_PROC 2)
 python_test(FILE dipolar_mdlc_p3m_scafacos_p2nfft.py MAX_NUM_PROC 1)
 python_test(FILE dipolar_direct_summation.py MAX_NUM_PROC 1 LABELS gpu)
 python_test(FILE dipolar_p3m.py MAX_NUM_PROC 2)