WIP: cast virtual sites into propagation scheme

src/core/EspressoSystemStandAlone.cpp

@@ -27,8 +27,6 @@
 #include "integrate.hpp"
 #include "magnetostatics/dipoles.hpp"
 #include "system/System.hpp"
-#include "virtual_sites.hpp"
-#include "virtual_sites/VirtualSitesOff.hpp"
 
 #include <utils/Vector.hpp>
 
@@ -46,9 +44,6 @@ EspressoSystemStandAlone::EspressoSystemStandAlone(int argc, char **argv) {
   Communication::init(mpi_env);
 
   // default-construct global state of the system
-#ifdef VIRTUAL_SITES
-  set_virtual_sites(std::make_shared<VirtualSitesOff>());
-#endif
   ::System::set_system(std::make_shared<::System::System>());
 }
 

src/core/event.cpp

@@ -48,7 +48,7 @@
 #include "particle_node.hpp"
 #include "system/System.hpp"
 #include "thermostat.hpp"
-#include "virtual_sites.hpp"
+#include "virtual_sites/relative.hpp"
 
 #include <utils/mpi/all_compare.hpp>
 
@@ -350,7 +350,7 @@ unsigned global_ghost_flags() {
 
 void update_dependent_particles() {
 #ifdef VIRTUAL_SITES
-  virtual_sites()->update();
+  vs_relative_update_particles(cell_structure);
   cells_update_ghosts(global_ghost_flags());
 #endif
 

src/core/forces.cpp

@@ -49,7 +49,7 @@
 #include "system/System.hpp"
 #include "thermostat.hpp"
 #include "thermostats/langevin_inline.hpp"
-#include "virtual_sites.hpp"
+#include "virtual_sites/relative.hpp"
 
 #include <boost/variant.hpp>
 
@@ -255,7 +255,7 @@ void force_calc(CellStructure &cell_structure, double time_step, double kT) {
 
 // VIRTUAL_SITES distribute forces
 #ifdef VIRTUAL_SITES
-  virtual_sites()->back_transfer_forces_and_torques();
+  vs_relative_back_transfer_forces_and_torques(cell_structure);
 #endif
 
   // Communication Step: ghost forces

src/core/integrate.cpp

@@ -55,7 +55,8 @@
 #include "rotation.hpp"
 #include "signalhandling.hpp"
 #include "thermostat.hpp"
-#include "virtual_sites.hpp"
+#include "virtual_sites/lb_tracers.hpp"
+#include "virtual_sites/relative.hpp"
 
 #include <boost/mpi/collectives/all_reduce.hpp>
 #include <boost/mpi/collectives/reduce.hpp>
@@ -340,7 +341,7 @@ int integrate(int n_steps, int reuse_forces) {
     lb_lbcoupling_deactivate();
 
 #ifdef VIRTUAL_SITES
-    virtual_sites()->update();
+    vs_relative_update_particles(cell_structure);
 #endif
 
     // Communication step: distribute ghost positions
@@ -418,7 +419,7 @@ int integrate(int n_steps, int reuse_forces) {
 #endif
 
 #ifdef VIRTUAL_SITES
-    virtual_sites()->update();
+    vs_relative_update_particles(cell_structure);
 #endif
 
     if (cell_structure.get_resort_particles() >= Cells::RESORT_LOCAL)
@@ -432,7 +433,7 @@ int integrate(int n_steps, int reuse_forces) {
     force_calc(cell_structure, time_step, temperature);
 
 #ifdef VIRTUAL_SITES
-    virtual_sites()->after_force_calc(time_step);
+    lb_tracers_add_particle_force_to_fluid(cell_structure, time_step);
 #endif
     integrator_step_2(particles, temperature);
     LeesEdwards::run_kernel<LeesEdwards::UpdateOffset>();
@@ -490,7 +491,7 @@ int integrate(int n_steps, int reuse_forces) {
       }
 
 #ifdef VIRTUAL_SITES
-      virtual_sites()->after_lb_propagation(time_step);
+      lb_tracers_propagate(cell_structure, time_step);
 #endif
 
 #ifdef COLLISION_DETECTION
@@ -523,7 +524,7 @@ int integrate(int n_steps, int reuse_forces) {
 #endif
 
 #ifdef VIRTUAL_SITES
-  virtual_sites()->update();
+  vs_relative_update_particles(cell_structure);
 #endif
 
   // Verlet list statistics

src/core/pressure.cpp

@@ -38,7 +38,7 @@
 #include "pressure_inline.hpp"
 #include "short_range_loop.hpp"
 #include "system/System.hpp"
-#include "virtual_sites.hpp"
+#include "virtual_sites/relative.hpp"
 
 #include "config/config.hpp"
 
@@ -111,7 +111,8 @@ std::shared_ptr<Observable_stat> calculate_pressure() {
 
 #ifdef VIRTUAL_SITES
   if (!obs_pressure.virtual_sites.empty()) {
-    auto const vs_pressure = virtual_sites()->pressure_tensor();
+    auto const vs_pressure =
+        vs_relative_pressure_tensor(cell_structure.local_particles());
     boost::copy(flatten(vs_pressure), obs_pressure.virtual_sites.begin());
   }
 #endif

src/core/virtual_sites.cpp

@@ -40,32 +40,20 @@
 #include <cstdio>
 #include <tuple>
 
-namespace {
-std::shared_ptr<VirtualSites> m_virtual_sites;
-}
-
-const std::shared_ptr<VirtualSites> &virtual_sites() { return m_virtual_sites; }
-
-void set_virtual_sites(std::shared_ptr<VirtualSites> const &v) {
-  m_virtual_sites = v;
-  recalc_forces = true;
-}
-
 #ifdef VIRTUAL_SITES_RELATIVE
 
 /** Calculate the rotation quaternion and distance between two particles */
 std::tuple<Utils::Quaternion<double>, double>
-calculate_vs_relate_to_params(Particle const &p_vs,
-                              Particle const &p_relate_to) {
+calculate_vs_relate_to_params(Particle const &p_vs, Particle const &p_relate_to,
+                              bool override_cutoff_check) {
   // get the distance between the particles
   auto d = ::box_geo.get_mi_vector(p_vs.pos(), p_relate_to.pos());
 
   // Check if the distance between virtual and non-virtual particles is larger
   // than minimum global cutoff. If so, warn user.
   auto const dist = d.norm();
   auto const min_global_cut = get_min_global_cut();
-  if (dist > min_global_cut && n_nodes > 1 &&
-      not virtual_sites()->get_override_cutoff_check()) {
+  if (dist > min_global_cut && n_nodes > 1 && not override_cutoff_check) {
     runtimeErrorMsg()
         << "Warning: The distance between virtual and non-virtual particle ("
         << dist << ") is larger than the minimum global cutoff ("

src/core/virtual_sites.hpp

@@ -24,16 +24,6 @@
 #ifdef VIRTUAL_SITES
 #include "Particle.hpp"
 
-#include "virtual_sites/VirtualSites.hpp"
-
-#include <memory>
-
-/** @brief Get active virtual sites implementation */
-const std::shared_ptr<VirtualSites> &virtual_sites();
-
-/** @brief Set active virtual sites implementation */
-void set_virtual_sites(std::shared_ptr<VirtualSites> const &v);
-
 #ifdef VIRTUAL_SITES_RELATIVE
 
 #include <utils/quaternion.hpp>
@@ -42,7 +32,8 @@ void set_virtual_sites(std::shared_ptr<VirtualSites> const &v);
 
 std::tuple<Utils::Quaternion<double>, double>
 calculate_vs_relate_to_params(Particle const &p_current,
-                              Particle const &p_relate_to);
+                              Particle const &p_relate_to,
+                              bool override_cutoff_check = false);
 
 /**
  * @brief Setup a virtual site to track a real particle.

src/core/virtual_sites/CMakeLists.txt

@@ -17,7 +17,5 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 
-target_sources(
-  espresso_core
-  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/VirtualSitesInertialessTracers.cpp
-          ${CMAKE_CURRENT_SOURCE_DIR}/VirtualSitesRelative.cpp)
+target_sources(espresso_core PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/lb_tracers.cpp
+                                     ${CMAKE_CURRENT_SOURCE_DIR}/relative.cpp)

src/core/virtual_sites/VirtualSitesInertialessTracers.cpp → src/core/virtual_sites/lb_tracers.cpp

@@ -20,8 +20,7 @@
 
 #ifdef VIRTUAL_SITES_INERTIALESS_TRACERS
 
-#include "VirtualSitesInertialessTracers.hpp"
-
+#include "PropagationModes.hpp"
 #include "cells.hpp"
 #include "errorhandling.hpp"
 #include "forces.hpp"
@@ -40,22 +39,27 @@ static bool lb_active_check() {
   return true;
 }
 
-void VirtualSitesInertialessTracers::after_force_calc(double time_step) {
+static bool is_lb_tracer(const Particle &p) {
+  return (p.propagation() & PropagationMode::TRANS_LB_TRACER);
+}
+
+void lb_tracers_add_particle_force_to_fluid(CellStructure &cell_struct,
+                                            double time_step) {
   auto const to_lb_units =
       (lattice_switch == ActiveLB::NONE) ? 0. : 1. / LB::get_agrid();
 
   // Distribute summed-up forces from physical particles to ghosts
-  init_forces_ghosts(cell_structure.ghost_particles());
+  init_forces_ghosts(cell_struct.ghost_particles());
   cells_update_ghosts(Cells::DATA_PART_FORCE);
 
   // Set to store ghost particles (ids) that have already been coupled
   LB::CouplingBookkeeping bookkeeping{};
   // Apply particle forces to the LB fluid at particle positions
   // For physical particles, also set particle velocity = fluid velocity
   for (auto const &particle_range :
-       {cell_structure.local_particles(), cell_structure.ghost_particles()}) {
+       {cell_struct.local_particles(), cell_struct.ghost_particles()}) {
     for (auto const &p : particle_range) {
-      if (!p.is_virtual())
+      if (!is_lb_tracer(p))
         continue;
       if (!lb_active_check()) {
         return;
@@ -69,16 +73,16 @@ void VirtualSitesInertialessTracers::after_force_calc(double time_step) {
   }
 
   // Clear ghost forces to avoid double counting later
-  init_forces_ghosts(cell_structure.ghost_particles());
+  init_forces_ghosts(cell_struct.ghost_particles());
 }
 
-void VirtualSitesInertialessTracers::after_lb_propagation(double time_step) {
+void lb_tracers_propagate(CellStructure &cell_struct, double time_step) {
   auto const to_md_units =
       (lattice_switch == ActiveLB::NONE) ? 0. : LB::get_lattice_speed();
 
   // Advect particles
-  for (auto &p : cell_structure.local_particles()) {
-    if (!p.is_virtual())
+  for (auto &p : cell_struct.local_particles()) {
+    if (!is_lb_tracer(p))
       continue;
     if (!lb_active_check()) {
       return;
@@ -91,7 +95,7 @@ void VirtualSitesInertialessTracers::after_lb_propagation(double time_step) {
     }
     // Verlet list update check
     if ((p.pos() - p.pos_at_last_verlet_update()).norm2() > skin * skin) {
-      cell_structure.set_resort_particles(Cells::RESORT_LOCAL);
+      cell_struct.set_resort_particles(Cells::RESORT_LOCAL);
     }
   }
 }