checkpointing

src/core/forces.cpp

@@ -211,19 +211,7 @@ void System::System::calculate_forces() {
                         collision_detection_cutoff});
 
   constraints->add_forces(particles, get_sim_time());
-
-  for (int i = 0; i < max_oif_objects; i++) {
-    // There are two global quantities that need to be evaluated:
-    // object's surface and object's volume.
-    auto const area_volume = boost::mpi::all_reduce(
-        comm_cart, calc_oif_global(i, *box_geo, *cell_structure), std::plus());
-    auto const oif_part_area = std::abs(area_volume[0]);
-    auto const oif_part_vol = std::abs(area_volume[1]);
-    if (oif_part_area < 1e-100 and oif_part_vol < 1e-100) {
-      break;
-    }
-    add_oif_global_forces(area_volume, i, *box_geo, *cell_structure);
-  }
+  oif_global->calculate_forces();
 
   // Must be done here. Forces need to be ghost-communicated
   immersed_boundaries->volume_conservation(*cell_structure);

src/core/immersed_boundary/ibm_tribend.cpp

@@ -95,6 +95,9 @@ IBMTribend::calc_forces(BoxGeometry const &box_geo, Particle const &p1,
 
 void IBMTribend::initialize(BoxGeometry const &box_geo,
                             CellStructure const &cell_structure) {
+  if (is_initialized) {
+    return;
+  }
   // Compute theta0
   if (flat) {
     theta0 = 0.;
@@ -128,4 +131,5 @@ void IBMTribend::initialize(BoxGeometry const &box_geo,
     if (desc < 0.)
       theta0 = 2. * std::numbers::pi - theta0;
   }
+  is_initialized = true;
 }

src/core/immersed_boundary/ibm_tribend.hpp

@@ -48,6 +48,7 @@ struct IBMTribend {
   std::tuple<int, int, int, int> p_ids;
 
   bool flat;
+  bool is_initialized;
 
   double cutoff() const { return 0.; }
 
@@ -61,7 +62,8 @@ struct IBMTribend {
                   CellStructure const &cell_structure);
 
   IBMTribend(int ind1, int ind2, int ind3, int ind4, double kb, bool flat)
-      : kb{kb}, theta0{0.}, p_ids{ind1, ind2, ind3, ind4}, flat{flat} {}
+      : kb{kb}, theta0{0.}, p_ids{ind1, ind2, ind3, ind4}, flat{flat},
+        is_initialized{false} {}
 
   /** Calculate the forces
    *  The equations can be found in Appendix C of @cite kruger12a.

src/core/immersed_boundary/ibm_triel.cpp

@@ -193,6 +193,9 @@ IBMTriel::calc_forces(Utils::Vector3d const &vec1,
 
 void IBMTriel::initialize(BoxGeometry const &box_geo,
                           CellStructure const &cell_structure) {
+  if (is_initialized) {
+    return;
+  }
   // collect particles from nodes
   auto const [ind1, ind2, ind3] = p_ids;
   auto const pos1 = get_ibm_particle_position(cell_structure, ind1);
@@ -221,4 +224,5 @@ void IBMTriel::initialize(BoxGeometry const &box_geo,
   b1 = (l0 * cosPhi0 - lp0) / area2;
   b2 = -(l0 * cosPhi0) / area2;
   area0 = 0.5 * area2;
+  is_initialized = true;
 }

src/core/immersed_boundary/ibm_triel.hpp

@@ -56,6 +56,7 @@ struct IBMTriel {
 
   /** Particle ids */
   std::tuple<int, int, int> p_ids;
+  bool is_initialized;
 
   double cutoff() const { return maxDist; }
 
@@ -71,7 +72,7 @@ struct IBMTriel {
            double k1, double k2)
       : l0{0.}, lp0{0.}, sinPhi0{0.}, cosPhi0{0.}, area0{0.}, a1{0.}, a2{0.},
         b1{0.}, b2{0.}, maxDist{maxDist}, elasticLaw{elasticLaw}, k1{k1},
-        k2{k2}, p_ids{ind1, ind2, ind3} {}
+        k2{k2}, p_ids{ind1, ind2, ind3}, is_initialized{false} {}
 
   /** Calculate the forces.
    *  The equations can be found in Appendix C of @cite kruger12a.

src/core/object-in-fluid/oif_global_forces.cpp

@@ -21,68 +21,68 @@
 
 #include "BoxGeometry.hpp"
 #include "Particle.hpp"
+#include "bonded_interactions/bonded_interaction_data.hpp"
 #include "cell_system/CellStructure.hpp"
+#include "communication.hpp"
+#include "oif_global_forces_params.hpp"
 #include "system/System.hpp"
 
-#include "bonded_interactions/bonded_interaction_data.hpp"
-
 #include <utils/Vector.hpp>
 #include <utils/math/triangle_functions.hpp>
 
+#include <boost/mpi/collectives/all_reduce.hpp>
+#include <boost/serialization/utility.hpp>
+
+#include <cassert>
+#include <cmath>
+#include <functional>
 #include <span>
 
-int max_oif_objects = 0;
+/** Calculate the mesh volume and area. */
+static auto calc_oif_mesh(int molType, BoxGeometry const &box_geo,
+                          CellStructure &cs,
+                          BondedInteractionsMap const &bonded_ias) {
 
-Utils::Vector2d calc_oif_global(int molType, BoxGeometry const &box_geo,
-                                CellStructure &cs) {
-  auto const &bonded_ias = *::System::get_system().bonded_ias;
-  // first-fold-then-the-same approach
-  double partArea = 0.0;
-  // z volume
-  double VOL_partVol = 0.;
+  double area = 0.;
+  double volume = 0.;
 
-  cs.bond_loop([&partArea, &VOL_partVol, &box_geo, &bonded_ias, molType](
+  cs.bond_loop([&area, &volume, &box_geo, &bonded_ias, molType](
                    Particle &p1, int bond_id, std::span<Particle *> partners) {
     if (p1.mol_id() != molType)
       return false;
 
-    if (boost::get<OifGlobalForcesBond>(bonded_ias.at(bond_id).get()) !=
-        nullptr) {
-      // remaining neighbors fetched
+    if (boost::get<OifGlobalForcesBond>(bonded_ias.at(bond_id).get())) {
       auto const p11 = box_geo.unfolded_position(p1.pos(), p1.image_box());
       auto const p22 = p11 + box_geo.get_mi_vector(partners[0]->pos(), p11);
       auto const p33 = p11 + box_geo.get_mi_vector(partners[1]->pos(), p11);
 
-      // unfolded positions correct
       auto const VOL_A = Utils::area_triangle(p11, p22, p33);
-      partArea += VOL_A;
+      area += VOL_A;
 
       auto const VOL_norm = Utils::get_n_triangle(p11, p22, p33);
       auto const VOL_dn = VOL_norm.norm();
-      auto const VOL_hz = 1.0 / 3.0 * (p11[2] + p22[2] + p33[2]);
-      VOL_partVol -= VOL_A * VOL_norm[2] / VOL_dn * VOL_hz;
+      auto const VOL_hz = (1.0 / 3.0) * (p11[2] + p22[2] + p33[2]);
+      volume -= VOL_A * VOL_norm[2] / VOL_dn * VOL_hz;
     }
 
     return false;
   });
 
-  return {{partArea, VOL_partVol}};
+  return Utils::Vector2d{{area, volume}};
 }
 
-void add_oif_global_forces(Utils::Vector2d const &area_volume, int molType,
-                           BoxGeometry const &box_geo, CellStructure &cs) {
-  auto const &bonded_ias = *::System::get_system().bonded_ias;
-  // first-fold-then-the-same approach
-  double area = area_volume[0];
-  double VOL_volume = area_volume[1];
+/** Distribute the OIF global forces to all particles in the mesh. */
+static void add_oif_global_forces(double area, double volume, int molType,
+                                  BoxGeometry const &box_geo, CellStructure &cs,
+                                  BondedInteractionsMap const &bonded_ias) {
 
-  cs.bond_loop([&box_geo, &bonded_ias, area, VOL_volume, molType](
+  cs.bond_loop([&box_geo, &bonded_ias, area, volume, molType](
                    Particle &p1, int bond_id, std::span<Particle *> partners) {
     if (p1.mol_id() != molType)
       return false;
 
-    if (auto const *iaparams =
-            boost::get<OifGlobalForcesBond>(bonded_ias.at(bond_id).get())) {
+    auto const *bond_ptr = bonded_ias.at(bond_id).get();
+    if (auto const *bond = boost::get<OifGlobalForcesBond>(bond_ptr)) {
       auto const p11 = box_geo.unfolded_position(p1.pos(), p1.image_box());
       auto const p22 = p11 + box_geo.get_mi_vector(partners[0]->pos(), p11);
       auto const p33 = p11 + box_geo.get_mi_vector(partners[1]->pos(), p11);
@@ -91,14 +91,13 @@ void add_oif_global_forces(Utils::Vector2d const &area_volume, int molType,
       // starting code from volume force
       auto const VOL_norm = Utils::get_n_triangle(p11, p22, p33).normalize();
       auto const VOL_A = Utils::area_triangle(p11, p22, p33);
-      auto const VOL_vv = (VOL_volume - iaparams->V0) / iaparams->V0;
+      auto const VOL_vv = (volume - bond->V0) / bond->V0;
 
-      auto const VOL_force =
-          (1.0 / 3.0) * iaparams->kv * VOL_vv * VOL_A * VOL_norm;
+      auto const VOL_force = (1. / 3.) * bond->kv * VOL_vv * VOL_A * VOL_norm;
 
       auto const h = (1. / 3.) * (p11 + p22 + p33);
 
-      auto const deltaA = (area - iaparams->A0_g) / iaparams->A0_g;
+      auto const deltaA = (area - bond->A0_g) / bond->A0_g;
 
       auto const m1 = h - p11;
       auto const m2 = h - p22;
@@ -108,7 +107,7 @@ void add_oif_global_forces(Utils::Vector2d const &area_volume, int molType,
       auto const m2_length = m2.norm();
       auto const m3_length = m3.norm();
 
-      auto const fac = iaparams->ka_g * VOL_A * deltaA /
+      auto const fac = bond->ka_g * VOL_A * deltaA /
                        (m1_length * m1_length + m2_length * m2_length +
                         m3_length * m3_length);
 
@@ -120,3 +119,22 @@ void add_oif_global_forces(Utils::Vector2d const &area_volume, int molType,
     return false;
   });
 }
+
+void OifGlobal::run_force_kernel() {
+  auto &system = get_system();
+  auto &box_geo = *system.box_geo;
+  auto &bonded_ias = *system.bonded_ias;
+  auto &cell_structure = *system.cell_structure;
+  for (int i = 0; i < max_oif_objects; ++i) {
+    // There are two global quantities that need to be evaluated:
+    // object's surface and object's volume.
+    auto const local = calc_oif_mesh(i, box_geo, cell_structure, bonded_ias);
+    auto const global = boost::mpi::all_reduce(comm_cart, local, std::plus());
+    auto const area = std::abs(global[0]);
+    auto const volume = std::abs(global[1]);
+    if (area < 1e-100 and volume < 1e-100) {
+      break;
+    }
+    add_oif_global_forces(area, volume, i, box_geo, cell_structure, bonded_ias);
+  }
+}

src/core/object-in-fluid/oif_global_forces.hpp

@@ -16,33 +16,26 @@
  * 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 CORE_OBJECT_IN_FLUID_OIF_GLOBAL_FORCES_HPP
-#define CORE_OBJECT_IN_FLUID_OIF_GLOBAL_FORCES_HPP
+
+#pragma once
+
 /** \file
- *  Routines to calculate the OIF global forces energy or/and and force
+ *  Routines to calculate the OIF global forces
  *  for a particle triple (triangle from mesh). See @cite dupin07a.
  */
 
-#include "BoxGeometry.hpp"
-#include "cell_system/CellStructure.hpp"
-#include "oif_global_forces_params.hpp"
-
-#include <utils/Vector.hpp>
-
-/** Calculate the OIF global force.
- *  Called in force_calc() from within forces.cpp
- *  - calculates the global area and global volume for a cell before the forces
- *    are handled
- *  - MPI synchronization with all reduce
- *  - !!! loop over particles from regular_decomposition !!!
- */
-Utils::Vector2d calc_oif_global(int molType, BoxGeometry const &box_geo,
-                                CellStructure &cs);
+#include "system/Leaf.hpp"
 
-/** Distribute the OIF global forces to all particles in the mesh. */
-void add_oif_global_forces(Utils::Vector2d const &area_volume, int molType,
-                           BoxGeometry const &box_geo, CellStructure &cs);
+class OifGlobal : public System::Leaf<OifGlobal> {
+public:
+  int max_oif_objects = 0;
 
-extern int max_oif_objects;
+  void calculate_forces() {
+    if (max_oif_objects > 0) {
+      run_force_kernel();
+    }
+  }
 
-#endif
+private:
+  void run_force_kernel();
+};

src/core/object-in-fluid/oif_global_forces_params.hpp

@@ -16,10 +16,8 @@
  * 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 OIF_GLOBAL_FORCES_PARAMS_HPP
-#define OIF_GLOBAL_FORCES_PARAMS_HPP
 
-#include <utils/Vector.hpp>
+#pragma once
 
 /** Parameters for OIF global forces
  *
@@ -46,16 +44,4 @@ struct OifGlobalForcesBond {
     this->V0 = V0;
     this->kv = kv;
   }
-
-private:
-  friend boost::serialization::access;
-  template <typename Archive>
-  void serialize(Archive &ar, long int /* version */) {
-    ar & A0_g;
-    ar & ka_g;
-    ar & V0;
-    ar & kv;
-  }
 };
-
-#endif

src/core/object-in-fluid/oif_local_forces.hpp

@@ -16,8 +16,8 @@
  * 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 CORE_OBJECT_IN_FLUID_OIF_LOCAL_FORCES_HPP
-#define CORE_OBJECT_IN_FLUID_OIF_LOCAL_FORCES_HPP
+
+#pragma once
 
 /** \file
  *  Routines to calculate the OIF local forces for a particle quadruple
@@ -80,21 +80,6 @@ struct OifLocalForcesBond {
   std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
   calc_forces(BoxGeometry const &box_geo, Particle const &p2,
               Particle const &p1, Particle const &p3, Particle const &p4) const;
-
-private:
-  friend boost::serialization::access;
-  template <typename Archive>
-  void serialize(Archive &ar, long int /* version */) {
-    ar & r0;
-    ar & ks;
-    ar & kslin;
-    ar & phi0;
-    ar & kb;
-    ar & A01;
-    ar & A02;
-    ar & kal;
-    ar & kvisc;
-  }
 };
 
 /** Compute the OIF local forces.
@@ -252,5 +237,3 @@ OifLocalForcesBond::calc_forces(BoxGeometry const &box_geo, Particle const &p2,
   }
   return std::make_tuple(force2, force1, force3, force4);
 }
-
-#endif

src/core/system/System.cpp

@@ -72,6 +72,7 @@ System::System(Private) {
   nonbonded_ias = std::make_shared<InteractionsNonBonded>();
   comfixed = std::make_shared<ComFixed>();
   galilei = std::make_shared<Galilei>();
+  oif_global = std::make_shared<OifGlobal>();
   immersed_boundaries = std::make_shared<ImmersedBoundaries>();
 #ifdef COLLISION_DETECTION
   collision_detection = std::make_shared<CollisionDetection>();
@@ -96,6 +97,8 @@ void System::initialize() {
   bonded_ias->bind_system(handle);
   thermostat->bind_system(handle);
   nonbonded_ias->bind_system(handle);
+  oif_global->bind_system(handle);
+  immersed_boundaries->bind_system(handle);
 #ifdef COLLISION_DETECTION
   collision_detection->bind_system(handle);
 #endif

src/core/system/System.hpp

@@ -48,6 +48,7 @@ class Thermostat;
 class ComFixed;
 class Galilei;
 class Observable_stat;
+class OifGlobal;
 class ImmersedBoundaries;
 class CollisionDetection;
 namespace BondBreakage {
@@ -280,6 +281,7 @@ class System : public std::enable_shared_from_this<System> {
   std::shared_ptr<Thermostat::Thermostat> thermostat;
   std::shared_ptr<ComFixed> comfixed;
   std::shared_ptr<Galilei> galilei;
+  std::shared_ptr<OifGlobal> oif_global;
   std::shared_ptr<ImmersedBoundaries> immersed_boundaries;
 #ifdef COLLISION_DETECTION
   std::shared_ptr<CollisionDetection> collision_detection;

src/core/system/System.impl.hpp

@@ -37,6 +37,7 @@
 #include "integrators/Propagation.hpp"
 #include "lees_edwards/lees_edwards.hpp"
 #include "nonbonded_interactions/nonbonded_interaction_data.hpp"
+#include "object-in-fluid/oif_global_forces.hpp"
 #include "thermostat.hpp"
 
 #include "BoxGeometry.hpp"