core: Encapsulate magnetostatics globals

The magnetostatics solver is now a member of the System class.
espressomd · Jul 5, 2023 · e4d68bc · e4d68bc
1 parent 03f0ac0
commit e4d68bc
Show file tree

Hide file tree

Showing 13 changed files with 216 additions and 142 deletions.
diff --git a/src/core/energy.cpp b/src/core/energy.cpp
@@ -66,7 +66,7 @@ std::shared_ptr<Observable_stat> calculate_energy() {
 
   auto const &system = System::get_system();
   auto const coulomb_kernel = system.coulomb.pair_energy_kernel();
-  auto const dipoles_kernel = Dipoles::pair_energy_kernel();
+  auto const dipoles_kernel = system.dipoles.pair_energy_kernel();
 
   short_range_loop(
       [&obs_energy, coulomb_kernel_ptr = coulomb_kernel.get_ptr()](
@@ -96,7 +96,7 @@ std::shared_ptr<Observable_stat> calculate_energy() {
 
 #ifdef DIPOLES
   /* calculate k-space part of magnetostatic interaction. */
-  obs_energy.dipolar[1] = Dipoles::calc_energy_long_range(local_parts);
+  obs_energy.dipolar[1] = system.dipoles.calc_energy_long_range(local_parts);
 #endif
 
   Constraints::constraints.add_energy(local_parts, get_sim_time(), obs_energy);
@@ -156,7 +156,8 @@ double particle_short_range_energy_contribution(int pid) {
 #ifdef DIPOLE_FIELD_TRACKING
 void calc_long_range_fields() {
   auto particles = cell_structure.local_particles();
-  Dipoles::calc_long_range_field(particles);
+  auto const &dipoles = System::get_system().dipoles;
+  dipoles.calc_long_range_field(particles);
 }
 
 REGISTER_CALLBACK(calc_long_range_fields)

diff --git a/src/core/event.cpp b/src/core/event.cpp
@@ -56,10 +56,6 @@
 
 /** whether the thermostat has to be reinitialized before integration */
 static bool reinit_thermo = true;
-#ifdef DIPOLES
-/** whether magnetostatics actor has to be reinitialized on observable calc */
-static bool reinit_magnetostatics = false;
-#endif
 
 void on_program_start() {
 #ifdef CUDA
@@ -125,7 +121,7 @@ void on_integration_start(double time_step) {
 #endif
 #ifdef DIPOLES
   {
-    auto const &actor = magnetostatics_actor;
+    auto const &actor = System::get_system().dipoles.solver;
     if (!Utils::Mpi::all_compare(comm_cart, static_cast<bool>(actor)) or
         (actor and !Utils::Mpi::all_compare(comm_cart, (*actor).which())))
       runtimeErrorMsg() << "Nodes disagree about dipolar long-range method";
@@ -147,11 +143,8 @@ void on_observable_calc() {
 #endif
 
 #ifdef DIPOLES
-  if (reinit_magnetostatics) {
-    Dipoles::on_observable_calc();
-    reinit_magnetostatics = false;
-  }
-#endif /* DIPOLES */
+  System::get_system().dipoles.on_observable_calc();
+#endif
 
   clear_particle_node();
 }
@@ -178,7 +171,9 @@ void on_particle_change() {
   }
 #endif
 #ifdef DIPOLES
-  reinit_magnetostatics = true;
+  if (System::is_system_set()) {
+    System::get_system().dipoles.on_particle_change();
+  }
 #endif
   recalc_forces = true;
 
@@ -194,8 +189,7 @@ void on_coulomb_and_dipoles_change() {
   System::get_system().coulomb.on_coulomb_change();
 #endif
 #ifdef DIPOLES
-  reinit_magnetostatics = true;
-  Dipoles::on_dipoles_change();
+  System::get_system().dipoles.on_dipoles_change();
 #endif
   on_short_range_ia_change();
 }
@@ -209,8 +203,7 @@ void on_coulomb_change() {
 
 void on_dipoles_change() {
 #ifdef DIPOLES
-  reinit_magnetostatics = true;
-  Dipoles::on_dipoles_change();
+  System::get_system().dipoles.on_dipoles_change();
 #endif
   on_short_range_ia_change();
 }
@@ -242,7 +235,7 @@ void on_boxl_change(bool skip_method_adaption) {
 #endif
 
 #ifdef DIPOLES
-    Dipoles::on_boxl_change();
+    System::get_system().dipoles.on_boxl_change();
 #endif
 
     LB::init();
@@ -268,7 +261,7 @@ void on_cell_structure_change() {
     System::get_system().coulomb.on_cell_structure_change();
 #endif
 #ifdef DIPOLES
-    Dipoles::on_cell_structure_change();
+    System::get_system().dipoles.on_cell_structure_change();
 #endif
   }
 #endif
@@ -286,7 +279,7 @@ void on_periodicity_change() {
 #endif
 
 #ifdef DIPOLES
-  Dipoles::on_periodicity_change();
+  System::get_system().dipoles.on_periodicity_change();
 #endif
 
 #ifdef STOKESIAN_DYNAMICS
@@ -321,7 +314,7 @@ void on_node_grid_change() {
   System::get_system().coulomb.on_node_grid_change();
 #endif
 #ifdef DIPOLES
-  Dipoles::on_node_grid_change();
+  System::get_system().dipoles.on_node_grid_change();
 #endif
   cells_re_init(cell_structure.decomposition_type());
 }

diff --git a/src/core/forces.cpp b/src/core/forces.cpp
@@ -170,7 +170,7 @@ void force_calc(CellStructure &cell_structure, double time_step, double kT) {
 
   auto const elc_kernel = espresso_system.coulomb.pair_force_elc_kernel();
   auto const coulomb_kernel = espresso_system.coulomb.pair_force_kernel();
-  auto const dipoles_kernel = Dipoles::pair_force_kernel();
+  auto const dipoles_kernel = espresso_system.dipoles.pair_force_kernel();
 
 #ifdef ELECTROSTATICS
   auto const coulomb_cutoff = espresso_system.coulomb.cutoff();
@@ -179,7 +179,7 @@ void force_calc(CellStructure &cell_structure, double time_step, double kT) {
 #endif
 
 #ifdef DIPOLES
-  auto const dipole_cutoff = Dipoles::cutoff();
+  auto const dipole_cutoff = espresso_system.dipoles.cutoff();
 #else
   auto const dipole_cutoff = INACTIVE_CUTOFF;
 #endif
@@ -261,7 +261,7 @@ void calc_long_range_forces(const ParticleRange &particles) {
 
 #ifdef DIPOLES
   /* calculate k-space part of the magnetostatic interaction. */
-  Dipoles::calc_long_range_force(particles);
+  Dipoles::get_dipoles().calc_long_range_force(particles);
 #endif // DIPOLES
 }
 

diff --git a/src/core/interactions.cpp b/src/core/interactions.cpp
@@ -42,7 +42,7 @@ static double recalc_long_range_cutoff() {
     max_cut_long_range = std::max(max_cut_long_range, system.coulomb.cutoff());
 #endif
 #ifdef DIPOLES
-    max_cut_long_range = std::max(max_cut_long_range, Dipoles::cutoff());
+    max_cut_long_range = std::max(max_cut_long_range, system.dipoles.cutoff());
 #endif
   }
 #endif
@@ -76,7 +76,7 @@ bool long_range_interactions_sanity_checks() {
     system.coulomb.sanity_checks();
 #endif
 #ifdef DIPOLES
-    Dipoles::sanity_checks();
+    system.dipoles.sanity_checks();
 #endif
   } catch (std::runtime_error const &err) {
     runtimeErrorMsg() << err.what();

diff --git a/src/core/magnetostatics/dipoles.cpp b/src/core/magnetostatics/dipoles.cpp
@@ -32,6 +32,7 @@
 #include "grid.hpp"
 #include "integrate.hpp"
 #include "npt.hpp"
+#include "system/System.hpp"
 
 #include <utils/Vector.hpp>
 #include <utils/constants.hpp>
@@ -43,67 +44,61 @@
 #include <cstdio>
 #include <stdexcept>
 
-boost::optional<MagnetostaticsActor> magnetostatics_actor;
-
 namespace Dipoles {
 
-void sanity_checks() {
-  if (magnetostatics_actor) {
-    boost::apply_visitor([](auto &actor) { actor->sanity_checks(); },
-                         *magnetostatics_actor);
+Solver const &get_dipoles() { return System::get_system().dipoles; }
+
+void Solver::sanity_checks() const {
+  if (solver) {
+    boost::apply_visitor([](auto &actor) { actor->sanity_checks(); }, *solver);
   }
 }
 
-void on_dipoles_change() {
-  visit_active_actor_try_catch([](auto &actor) { actor->init(); },
-                               magnetostatics_actor);
+void Solver::on_dipoles_change() {
+  reinit_on_observable_calc = true;
+  visit_active_actor_try_catch([](auto &actor) { actor->init(); }, solver);
 }
 
-void on_boxl_change() {
+void Solver::on_boxl_change() {
   visit_active_actor_try_catch([](auto &actor) { actor->on_boxl_change(); },
-                               magnetostatics_actor);
+                               solver);
 }
 
-void on_node_grid_change() {
-  if (magnetostatics_actor) {
+void Solver::on_node_grid_change() {
+  if (solver) {
     boost::apply_visitor([](auto &actor) { actor->on_node_grid_change(); },
-                         *magnetostatics_actor);
+                         *solver);
   }
 }
 
-void on_periodicity_change() {
+void Solver::on_periodicity_change() {
   visit_active_actor_try_catch(
-      [](auto &actor) { actor->on_periodicity_change(); },
-      magnetostatics_actor);
+      [](auto &actor) { actor->on_periodicity_change(); }, solver);
 }
 
-void on_cell_structure_change() {
+void Solver::on_cell_structure_change() {
   visit_active_actor_try_catch(
-      [](auto &actor) { actor->on_cell_structure_change(); },
-      magnetostatics_actor);
-}
-
-void calc_pressure_long_range() {
-  if (magnetostatics_actor) {
-    runtimeWarningMsg() << "pressure calculated, but pressure not implemented.";
-  }
+      [](auto &actor) { actor->on_cell_structure_change(); }, solver);
 }
 
-double cutoff() {
+double Solver::cutoff() const {
 #ifdef DP3M
-  if (auto dp3m = get_actor_by_type<DipolarP3M>(magnetostatics_actor)) {
+  if (auto dp3m = get_actor_by_type<DipolarP3M>(solver)) {
     return dp3m->dp3m.params.r_cut;
   }
 #endif
   return -1.;
 }
 
-void on_observable_calc() {
+void Solver::on_observable_calc() {
+  if (reinit_on_observable_calc) {
 #ifdef DP3M
-  if (auto dp3m = get_actor_by_type<DipolarP3M>(magnetostatics_actor)) {
-    dp3m->count_magnetic_particles();
-  }
+    if (auto dp3m = get_actor_by_type<DipolarP3M>(solver)) {
+      dp3m->count_magnetic_particles();
+    }
 #endif
+    reinit_on_observable_calc = false;
+  }
 }
 
 struct LongRangeForce : public boost::static_visitor<void> {
@@ -205,24 +200,29 @@ struct LongRangeField : public boost::static_visitor<void> {
 };
 #endif
 
-void calc_long_range_force(ParticleRange const &particles) {
-  if (magnetostatics_actor) {
-    boost::apply_visitor(LongRangeForce(particles), *magnetostatics_actor);
+void Solver::calc_pressure_long_range() const {
+  if (solver) {
+    runtimeWarningMsg() << "pressure calculated, but pressure not implemented.";
+  }
+}
+
+void Solver::calc_long_range_force(ParticleRange const &particles) const {
+  if (solver) {
+    boost::apply_visitor(LongRangeForce(particles), *solver);
   }
 }
 
-double calc_energy_long_range(ParticleRange const &particles) {
-  if (magnetostatics_actor) {
-    return boost::apply_visitor(LongRangeEnergy(particles),
-                                *magnetostatics_actor);
+double Solver::calc_energy_long_range(ParticleRange const &particles) const {
+  if (solver) {
+    return boost::apply_visitor(LongRangeEnergy(particles), *solver);
   }
   return 0.;
 }
 
 #ifdef DIPOLE_FIELD_TRACKING
-void calc_long_range_field(ParticleRange const &particles) {
-  if (magnetostatics_actor) {
-    boost::apply_visitor(LongRangeField(particles), *magnetostatics_actor);
+void Solver::calc_long_range_field(ParticleRange const &particles) const {
+  if (solver) {
+    boost::apply_visitor(LongRangeField(particles), *solver);
   }
 }
 #endif

diff --git a/src/core/magnetostatics/dipoles.hpp b/src/core/magnetostatics/dipoles.hpp
@@ -26,6 +26,8 @@
 
 #include "actor/traits.hpp"
 
+#include "magnetostatics/solver.hpp"
+
 #include "magnetostatics/barnes_hut_gpu.hpp"
 #include "magnetostatics/dipolar_direct_sum.hpp"
 #include "magnetostatics/dipolar_direct_sum_gpu.hpp"
@@ -47,24 +49,6 @@
 #include <stdexcept>
 #include <type_traits>
 
-using MagnetostaticsActor =
-    boost::variant<std::shared_ptr<DipolarDirectSum>,
-#ifdef DIPOLAR_DIRECT_SUM
-                   std::shared_ptr<DipolarDirectSumGpu>,
-#endif
-#ifdef DIPOLAR_BARNES_HUT
-                   std::shared_ptr<DipolarBarnesHutGpu>,
-#endif
-#ifdef DP3M
-                   std::shared_ptr<DipolarP3M>,
-#endif
-#ifdef SCAFACOS_DIPOLES
-                   std::shared_ptr<DipolarScafacos>,
-#endif
-                   std::shared_ptr<DipolarLayerCorrection>>;
-
-extern boost::optional<MagnetostaticsActor> magnetostatics_actor;
-
 /** Get the magnetostatics prefactor. */
 struct GetDipolesPrefactor : public boost::static_visitor<double> {
   template <typename T>
@@ -81,7 +65,6 @@ struct DipolesSanityChecks : public boost::static_visitor<void> {
 };
 
 namespace Dipoles {
-
 namespace traits {
 
 /** @brief Whether an actor is a solver. */
@@ -95,35 +78,7 @@ template <> struct has_dipole_fields<DipolarDirectSum> : std::true_type {};
 #endif // DIPOLE_FIELD_TRACKING
 
 } // namespace traits
-
-void calc_pressure_long_range();
-
-void sanity_checks();
-double cutoff();
-
-void on_observable_calc();
-void on_dipoles_change();
-void on_boxl_change();
-void on_node_grid_change();
-void on_periodicity_change();
-void on_cell_structure_change();
-
-void calc_long_range_force(ParticleRange const &particles);
-double calc_energy_long_range(ParticleRange const &particles);
-#ifdef DIPOLE_FIELD_TRACKING
-void calc_long_range_field(ParticleRange const &particles);
-#endif
-
-namespace detail {
-bool flag_all_reduce(bool flag);
-} // namespace detail
-
-} // namespace Dipoles
-#else // DIPOLES
-#include <cstddef>
-namespace Dipoles {
-constexpr std::size_t pressure_n() { return 0; }
-constexpr std::size_t energy_n() { return 0; }
 } // namespace Dipoles
+
 #endif // DIPOLES
 #endif