Merge #3053

3053: Max cut r=RudolfWeeber a=fweik

Description of changes:
 - Replaced some globals related to cutoff calculation
 - Store range limits in cell_structure
 - Consistently use `INACTIVE_CUTOFF`/`-1.` for cutoffs
   of disabled kernels.

Further refactoring is blocked by the fact that the short-range kernels
with different cutoffs are fused together, and by the unclear resposibility
for the verlet list. The performance hacks that mix everything together
also do not help.

Co-authored-by: Florian Weik <fweik@icp.uni-stuttgart.de>
Co-authored-by: Jean-Noël Grad <jgrad@icp.uni-stuttgart.de>

src/config/config.hpp

@@ -40,11 +40,6 @@
 #define ONEPART_DEBUG_ID 13
 #endif
 
-/** CELLS: Default value for the maximal number of cells per node. */
-#ifndef CELLS_MAX_NUM_CELLS
-#define CELLS_MAX_NUM_CELLS 32768
-#endif
-
 /** P3M: Default for number of interpolation points of the charge
     assignment function. */
 #ifndef P3M_N_INTERPOL

src/core/MpiCallbacks.hpp

@@ -131,6 +131,9 @@ template <class F, class... Args>
 struct callback_void_t final : public callback_concept_t {
   F m_f;
 
+  callback_void_t(callback_void_t const &) = delete;
+  callback_void_t(callback_void_t &&) = delete;
+
   template <class FRef>
   explicit callback_void_t(FRef &&f) : m_f(std::forward<FRef>(f)) {}
   void operator()(boost::mpi::communicator const &,
@@ -150,6 +153,9 @@ template <class F, class... Args>
 struct callback_one_rank_t final : public callback_concept_t {
   F m_f;
 
+  callback_one_rank_t(callback_one_rank_t const &) = delete;
+  callback_one_rank_t(callback_one_rank_t &&) = delete;
+
   template <class FRef>
   explicit callback_one_rank_t(FRef &&f) : m_f(std::forward<FRef>(f)) {}
   void operator()(boost::mpi::communicator const &comm,

src/core/bonded_interactions/bonded_interaction_data.cpp

@@ -1,46 +1,39 @@
 #include "bonded_interaction_data.hpp"
-#include "bonded_interactions/thermalized_bond.hpp"
 #include "communication.hpp"
 
 #include <utils/constants.hpp>
 
 std::vector<Bonded_ia_parameters> bonded_ia_params;
 
-void recalc_maximal_cutoff_bonded() {
-  int i;
-  double max_cut_tmp;
+double recalc_maximal_cutoff_bonded() {
+  auto max_cut_bonded = -1.;
 
-  max_cut_bonded = 0.0;
-
-  for (i = 0; i < bonded_ia_params.size(); i++) {
-    switch (bonded_ia_params[i].type) {
+  for (auto const &bonded_ia_param : bonded_ia_params) {
+    switch (bonded_ia_param.type) {
     case BONDED_IA_FENE:
-      max_cut_tmp =
-          bonded_ia_params[i].p.fene.r0 + bonded_ia_params[i].p.fene.drmax;
-      if (max_cut_bonded < max_cut_tmp)
-        max_cut_bonded = max_cut_tmp;
+      max_cut_bonded =
+          std::max(max_cut_bonded,
+                   bonded_ia_param.p.fene.r0 + bonded_ia_param.p.fene.drmax);
       break;
     case BONDED_IA_HARMONIC:
-      if ((bonded_ia_params[i].p.harmonic.r_cut > 0) &&
-          (max_cut_bonded < bonded_ia_params[i].p.harmonic.r_cut))
-        max_cut_bonded = bonded_ia_params[i].p.harmonic.r_cut;
+      max_cut_bonded =
+          std::max(max_cut_bonded, bonded_ia_param.p.harmonic.r_cut);
       break;
     case BONDED_IA_THERMALIZED_DIST:
-      if ((bonded_ia_params[i].p.thermalized_bond.r_cut > 0) &&
-          (max_cut_bonded < bonded_ia_params[i].p.thermalized_bond.r_cut))
-        max_cut_bonded = bonded_ia_params[i].p.thermalized_bond.r_cut;
+      max_cut_bonded =
+          std::max(max_cut_bonded, bonded_ia_param.p.thermalized_bond.r_cut);
       break;
     case BONDED_IA_RIGID_BOND:
-      if (max_cut_bonded < sqrt(bonded_ia_params[i].p.rigid_bond.d2))
-        max_cut_bonded = sqrt(bonded_ia_params[i].p.rigid_bond.d2);
+      max_cut_bonded =
+          std::max(max_cut_bonded, std::sqrt(bonded_ia_param.p.rigid_bond.d2));
       break;
     case BONDED_IA_TABULATED_DISTANCE:
-      if (max_cut_bonded < bonded_ia_params[i].p.tab.pot->cutoff())
-        max_cut_bonded = bonded_ia_params[i].p.tab.pot->cutoff();
+      max_cut_bonded =
+          std::max(max_cut_bonded, bonded_ia_param.p.tab.pot->cutoff());
       break;
     case BONDED_IA_IBM_TRIEL:
-      if (max_cut_bonded < bonded_ia_params[i].p.ibm_triel.maxDist)
-        max_cut_bonded = bonded_ia_params[i].p.ibm_triel.maxDist;
+      max_cut_bonded =
+          std::max(max_cut_bonded, bonded_ia_param.p.ibm_triel.maxDist);
       break;
     default:
       break;
@@ -49,19 +42,18 @@ void recalc_maximal_cutoff_bonded() {
 
   /* dihedrals: the central particle is indirectly connected to the fourth
    * particle via the third particle, so we have to double the cutoff */
-  max_cut_tmp = 2.0 * max_cut_bonded;
-  for (i = 0; i < bonded_ia_params.size(); i++) {
-    switch (bonded_ia_params[i].type) {
+  for (auto const &bonded_ia_param : bonded_ia_params) {
+    switch (bonded_ia_param.type) {
     case BONDED_IA_DIHEDRAL:
-      max_cut_bonded = max_cut_tmp;
-      break;
     case BONDED_IA_TABULATED_DIHEDRAL:
-      max_cut_bonded = max_cut_tmp;
+      max_cut_bonded *= 2;
       break;
     default:
       break;
     }
   }
+
+  return max_cut_bonded;
 }
 
 void make_bond_type_exist(int type) {

src/core/bonded_interactions/bonded_interaction_data.hpp

@@ -345,13 +345,6 @@ struct Bonded_ia_parameters {
 /** Field containing the parameters of the bonded ia types */
 extern std::vector<Bonded_ia_parameters> bonded_ia_params;
 
-/** @brief Maximal interaction cutoff for bonded interactions (real space).
- *  This value must be as large as the maximal interaction range in the list
- *  of bonded interactions. This is necessary to ensure that in a parallel
- *  simulation, a compute node has access to both bond partners.
- */
-extern double max_cut_bonded;
-
 /** Makes sure that \ref bonded_ia_params is large enough to cover the
  *  parameters for the bonded interaction type.
  *  Attention: 1: There is no initialization done here.
@@ -442,10 +435,8 @@ inline bool pair_bond_enum_exists_between(Particle const *const p1,
  *  is set to twice the maximal cutoff because the particle in which the bond
  *  is stored is only bonded to the first two partners, one of which has an
  *  additional bond to the third partner.
- *
- *  The result is stored in global variable @ref max_cut_bonded.
  */
-void recalc_maximal_cutoff_bonded();
+double recalc_maximal_cutoff_bonded();
 
 int virtual_set_params(int bond_type);
 #endif

src/core/bonded_interactions/bonded_interactions.dox

@@ -113,7 +113,7 @@
  *
  *  * In bonded_interaction_data.cpp:
  *    - in function @ref recalc_maximal_cutoff_bonded(): add a case for the new
- *      interaction which makes sure that @ref max_cut_bonded is as large as
+ *      interaction which makes sure that @ref max_cut is as large as
  *      the interaction range of the new interaction. This is only relevant to
  *      pairwise bonds and dihedral bonds. The range can be calculated by a
  *      formula, so it is not strictly necessary that the maximal interaction

src/core/cells.cpp

@@ -61,8 +61,6 @@ CellPList ghost_cells = {nullptr, 0, 0};
 /** Type of cell structure in use */
 CellStructure cell_structure;
 
-double max_range = 0.0;
-
 /** On of Cells::Resort, announces the level of resort needed.
  */
 unsigned resort_particles = Cells::RESORT_NONE;
@@ -186,24 +184,26 @@ static void topology_release(int cs) {
 }
 
 /** Choose the topology init function of a certain cell system. */
-void topology_init(int cs, CellPList *local) {
+void topology_init(int cs, double range, CellPList *local) {
   /** broadcast the flag for using Verlet list */
   boost::mpi::broadcast(comm_cart, cell_structure.use_verlet_list, 0);
 
   switch (cs) {
+  /* Default to DD */
   case CELL_STRUCTURE_NONEYET:
+    topology_init(CELL_STRUCTURE_DOMDEC, range, local);
     break;
   case CELL_STRUCTURE_CURRENT:
-    topology_init(cell_structure.type, local);
+    topology_init(cell_structure.type, range, local);
     break;
   case CELL_STRUCTURE_DOMDEC:
-    dd_topology_init(local, node_grid);
+    dd_topology_init(local, node_grid, range);
     break;
   case CELL_STRUCTURE_NSQUARE:
     nsq_topology_init(local);
     break;
   case CELL_STRUCTURE_LAYERED:
-    layered_topology_init(local, node_grid);
+    layered_topology_init(local, node_grid, range);
     break;
   default:
     fprintf(stderr,
@@ -248,7 +248,7 @@ static void invalidate_ghosts() {
 
 /************************************************************/
 
-void cells_re_init(int new_cs) {
+void cells_re_init(int new_cs, double range) {
   CellPList tmp_local;
 
   CELL_TRACE(fprintf(stderr, "%d: cells_re_init: convert type (%d->%d)\n",
@@ -264,7 +264,8 @@ void cells_re_init(int new_cs) {
   /* MOVE old cells to temporary buffer */
   auto tmp_cells = std::move(cells);
 
-  topology_init(new_cs, &tmp_local);
+  topology_init(new_cs, range, &tmp_local);
+  cell_structure.min_range = range;
 
   clear_particle_node();
 
@@ -275,8 +276,6 @@ void cells_re_init(int new_cs) {
     cell.resize(0);
   }
 
-  CELL_TRACE(fprintf(stderr, "%d: old cells deallocated\n", this_node));
-
   /* to enforce initialization of the ghost cells */
   resort_particles = Cells::RESORT_GLOBAL;
 
@@ -427,22 +426,17 @@ void cells_resort_particles(int global_flag) {
 /*************************************************/
 
 void cells_on_geometry_change(int flags) {
-  if (max_cut > 0.0) {
-    max_range = max_cut + skin;
-  } else
-    /* if no interactions yet, we also don't need a skin */
-    max_range = 0.0;
-
-  CELL_TRACE(fprintf(stderr, "%d: on_geometry_change with max range %f\n",
-                     this_node, max_range));
+  /* Consider skin only if there are actually interactions */
+  auto const range = (max_cut > 0.) ? max_cut + skin : INACTIVE_CUTOFF;
+  cell_structure.min_range = range;
 
   switch (cell_structure.type) {
   case CELL_STRUCTURE_DOMDEC:
-    dd_on_geometry_change(flags, node_grid);
+    dd_on_geometry_change(flags, node_grid, range);
     break;
   case CELL_STRUCTURE_LAYERED:
     /* there is no fast version, always redo everything. */
-    cells_re_init(CELL_STRUCTURE_LAYERED);
+    cells_re_init(CELL_STRUCTURE_LAYERED, range);
     break;
   case CELL_STRUCTURE_NSQUARE:
     break;

src/core/cells.hpp

@@ -147,6 +147,16 @@ struct CellStructure {
 
   bool use_verlet_list = true;
 
+  /** Maximal pair range supported by current
+   * cell system.
+   */
+  Utils::Vector3d max_range = {};
+
+  /**
+   * Minimum range that has to be supported.
+   */
+  double min_range;
+
   /** returns the global local_cells */
   CellPList local_cells() const;
   /** returns the global ghost_cells */
@@ -189,12 +199,6 @@ extern CellPList ghost_cells;
 /** Type of cell structure in use ( \ref Cell Structure ). */
 extern CellStructure cell_structure;
 
-/** Maximal interaction range - also the minimum cell size. Any
- *  cellsystem makes sure that the particle pair loop visits all pairs
- *  of particles that are closer than this.
- */
-extern double max_range;
-
 /** If non-zero, cell systems should reset the position for checking
  *  the Verlet criterion. Moreover, the Verlet list has to be
  *  rebuilt.
@@ -208,14 +212,12 @@ extern int rebuild_verletlist;
 /************************************************************/
 /*@{*/
 
-/** Switch for choosing the topology init function of a certain cell system. */
-void topology_init(int cs, CellPList *local);
-
 /** Reinitialize the cell structures.
  *  @param new_cs gives the new topology to use afterwards. May be set to
  *  \ref CELL_STRUCTURE_CURRENT for not changing it.
+ *  @param range Desired interaction range
  */
-void cells_re_init(int new_cs);
+void cells_re_init(int new_cs, double range);
 
 /** Reallocate the list of all cells (\ref cells::cells). */
 void realloc_cells(int size);

src/core/communication.cpp

@@ -73,6 +73,7 @@
 #include <utils/u32_to_u64.hpp>
 
 #include <boost/mpi.hpp>
+#include <boost/range/algorithm/min_element.hpp>
 #include <boost/serialization/array.hpp>
 #include <boost/serialization/string.hpp>
 #include <boost/serialization/utility.hpp>
@@ -394,7 +395,7 @@ void mpi_gather_stats(int job, void *result, void *result_t, void *result_nb,
   switch (job) {
   case 1:
     mpi_call(mpi_gather_stats_slave, -1, 1);
-    energy_calc((double *)result);
+    energy_calc((double *)result, sim_time);
     break;
   case 2:
     /* calculate and reduce (sum up) virials for 'analyze pressure' or
@@ -438,7 +439,7 @@ void mpi_gather_stats_slave(int, int job) {
   switch (job) {
   case 1:
     /* calculate and reduce (sum up) energies */
-    energy_calc(nullptr);
+    energy_calc(nullptr, sim_time);
     break;
   case 2:
     /* calculate and reduce (sum up) virials for 'analyze pressure' or 'analyze
@@ -549,10 +550,12 @@ void mpi_rescale_particles_slave(int, int dir) {
 
 void mpi_bcast_cell_structure(int cs) {
   mpi_call(mpi_bcast_cell_structure_slave, -1, cs);
-  cells_re_init(cs);
+  cells_re_init(cs, cell_structure.min_range);
 }
 
-void mpi_bcast_cell_structure_slave(int, int cs) { cells_re_init(cs); }
+void mpi_bcast_cell_structure_slave(int, int cs) {
+  cells_re_init(cs, cell_structure.min_range);
+}
 
 /*************** REQ_BCAST_NPTISO_GEOM *****************/