Introduce a System class in the core

src/core/EspressoSystemStandAlone.cpp

@@ -24,6 +24,7 @@
 #include "event.hpp"
 #include "grid.hpp"
 #include "integrate.hpp"
+#include "system/System.hpp"
 #include "virtual_sites.hpp"
 #include "virtual_sites/VirtualSitesOff.hpp"
 
@@ -46,6 +47,7 @@ EspressoSystemStandAlone::EspressoSystemStandAlone(int argc, char **argv) {
 #ifdef VIRTUAL_SITES
   set_virtual_sites(std::make_shared<VirtualSitesOff>());
 #endif
+  ::System::set_system(std::make_shared<::System::System>());
 }
 
 void EspressoSystemStandAlone::set_box_l(Utils::Vector3d const &box_l) const {

src/core/electrostatics/mmm1d_gpu.cpp

@@ -52,7 +52,7 @@ CoulombMMM1DGpu::CoulombMMM1DGpu(double prefactor, double maxPWerror,
     throw std::domain_error("Parameter 'bessel_cutoff' must be > 0");
   }
 
-  auto &gpu_particle_data = System::Instance().gpu;
+  auto &gpu_particle_data = System::get_system().gpu;
   gpu_particle_data.enable_property(GpuParticleData::prop::force);
   gpu_particle_data.enable_property(GpuParticleData::prop::pos);
   gpu_particle_data.enable_property(GpuParticleData::prop::q);
@@ -75,7 +75,7 @@ void CoulombMMM1DGpu::sanity_checks_cell_structure() const {
 }
 
 void CoulombMMM1DGpu::tune() {
-  System::Instance().gpu.init();
+  System::get_system().gpu.init();
   if (this_node == 0) {
     setup();
     tune(maxPWerror, far_switch_radius, bessel_cutoff);

src/core/electrostatics/mmm1d_gpu_cuda.cu

@@ -146,8 +146,8 @@ static auto numBlocks(std::size_t n_part) {
 }
 
 void CoulombMMM1DGpu::setup() {
-  auto &gpu_particle_data = System::Instance().gpu;
-  auto const box_z = static_cast<float>(System::Instance().box()[2]);
+  auto &gpu_particle_data = System::get_system().gpu;
+  auto const box_z = static_cast<float>(System::get_system().box()[2]);
   auto const n_part = gpu_particle_data.n_particles();
   if (not m_is_tuned and n_part != 0) {
     set_params(box_z, prefactor, maxPWerror, far_switch_radius, bessel_cutoff);
@@ -518,7 +518,7 @@ void CoulombMMM1DGpu::add_long_range_forces() {
     throw std::runtime_error("MMM1D was not initialized correctly");
   }
 
-  auto &gpu = System::Instance().gpu;
+  auto &gpu = System::get_system().gpu;
   auto const positions_device = gpu.get_particle_positions_device();
   auto const charges_device = gpu.get_particle_charges_device();
   auto const forces_device = gpu.get_particle_forces_device();
@@ -551,7 +551,7 @@ void CoulombMMM1DGpu::add_long_range_energy() {
     throw std::runtime_error("MMM1D was not initialized correctly");
   }
 
-  auto &gpu = System::Instance().gpu;
+  auto &gpu = System::get_system().gpu;
   auto const positions_device = gpu.get_particle_positions_device();
   auto const charges_device = gpu.get_particle_charges_device();
   auto const energy_device = gpu.get_energy_device();
@@ -573,7 +573,7 @@ float CoulombMMM1DGpu::force_benchmark() {
   float elapsedTime;
   float *dev_f_benchmark;
 
-  auto &gpu = System::Instance().gpu;
+  auto &gpu = System::get_system().gpu;
   auto const positions_device = gpu.get_particle_positions_device();
   auto const charges_device = gpu.get_particle_charges_device();
   auto const n_part = gpu.n_particles();

src/core/electrostatics/p3m_gpu.cpp

@@ -50,7 +50,7 @@ void CoulombP3MGPU::init() {
 void CoulombP3MGPU::init_cpu_kernels() { CoulombP3M::init(); }
 
 void CoulombP3MGPU::request_gpu() const {
-  auto &gpu_particle_data = System::Instance().gpu;
+  auto &gpu_particle_data = System::get_system().gpu;
   gpu_particle_data.enable_property(GpuParticleData::prop::force);
   gpu_particle_data.enable_property(GpuParticleData::prop::q);
   gpu_particle_data.enable_property(GpuParticleData::prop::pos);

src/core/electrostatics/p3m_gpu.hpp

@@ -39,6 +39,7 @@ struct CoulombP3MGPU : public CoulombP3M {
 
   void init();
   void on_activation() {
+    request_gpu();
     CoulombP3M::on_activation();
     if (is_tuned()) {
       init_cpu_kernels();

src/core/electrostatics/p3m_gpu_cuda.cu

@@ -552,8 +552,7 @@ void p3m_gpu_init(int cao, const int mesh[3], double alpha) {
   if (mesh[0] == -1 && mesh[1] == -1 && mesh[2] == -1)
     throw std::runtime_error("P3M: invalid mesh size");
 
-  auto &gpu_particle_data = System::Instance().gpu;
-
+  auto &gpu_particle_data = System::get_system().gpu;
   bool do_reinit = false, mesh_changed = false;
   p3m_gpu_data.n_part = static_cast<int>(gpu_particle_data.n_particles());
 
@@ -576,7 +575,7 @@ void p3m_gpu_init(int cao, const int mesh[3], double alpha) {
     do_reinit = true;
   }
 
-  auto const box_l = System::Instance().box();
+  auto const box_l = System::get_system().box();
 
   if (!p3m_gpu_data_initialized || (p3m_gpu_data.box[0] != box_l[0]) ||
       (p3m_gpu_data.box[1] != box_l[1]) || (p3m_gpu_data.box[2] != box_l[2])) {
@@ -680,7 +679,7 @@ void p3m_gpu_init(int cao, const int mesh[3], double alpha) {
  *  \brief The long-range part of the P3M algorithm.
  */
 void p3m_gpu_add_farfield_force(double prefactor) {
-  auto &gpu = System::Instance().gpu;
+  auto &gpu = System::get_system().gpu;
   p3m_gpu_data.n_part = static_cast<int>(gpu.n_particles());
 
   if (p3m_gpu_data.n_part == 0)

src/core/energy.cpp

@@ -51,7 +51,7 @@ std::shared_ptr<Observable_stat> calculate_energy() {
 
   auto &obs_energy = *obs_energy_ptr;
 #if defined(CUDA) and (defined(ELECTROSTATICS) or defined(DIPOLES))
-  auto &gpu_particle_data = System::Instance().gpu;
+  auto &gpu_particle_data = System::get_system().gpu;
   gpu_particle_data.clear_energy_on_device();
   gpu_particle_data.update();
 #endif

src/core/forces.cpp

@@ -146,7 +146,7 @@ void force_calc(CellStructure &cell_structure, double time_step, double kT) {
   ESPRESSO_PROFILER_CXX_MARK_FUNCTION;
 
 #ifdef CUDA
-  auto &espresso_system = System::Instance();
+  auto &espresso_system = System::get_system();
   espresso_system.gpu.update();
 #endif
 

src/core/magnetostatics/barnes_hut_gpu.cpp

@@ -42,7 +42,7 @@ DipolarBarnesHutGpu::DipolarBarnesHutGpu(double prefactor, double epssq,
   if (m_epssq <= 0.) {
     throw std::domain_error("Parameter 'epssq' must be > 0");
   }
-  auto &gpu_particle_data = System::Instance().gpu;
+  auto &gpu_particle_data = System::get_system().gpu;
   gpu_particle_data.enable_property(GpuParticleData::prop::force);
   gpu_particle_data.enable_property(GpuParticleData::prop::torque);
   gpu_particle_data.enable_property(GpuParticleData::prop::pos);
@@ -65,7 +65,7 @@ int call_kernel(void (*fp)(Args...), ArgRef &&...args) {
 }
 
 int DipolarBarnesHutGpu::initialize_data_structure() {
-  auto &gpu = System::Instance().gpu;
+  auto &gpu = System::get_system().gpu;
   auto const n_part = static_cast<int>(gpu.n_particles());
   auto const error_code = call_kernel(allocBHmemCopy, n_part, &m_bh_data);
 
@@ -84,7 +84,7 @@ int DipolarBarnesHutGpu::initialize_data_structure() {
 }
 
 void DipolarBarnesHutGpu::add_long_range_forces() {
-  auto &gpu_particle_data = System::Instance().gpu;
+  auto &gpu_particle_data = System::get_system().gpu;
   gpu_particle_data.update();
   if (this_node == 0) {
     if (initialize_data_structure() == ES_OK) {
@@ -97,7 +97,7 @@ void DipolarBarnesHutGpu::add_long_range_forces() {
 }
 
 void DipolarBarnesHutGpu::long_range_energy() {
-  auto &gpu_particle_data = System::Instance().gpu;
+  auto &gpu_particle_data = System::get_system().gpu;
   gpu_particle_data.update();
   if (this_node == 0) {
     if (initialize_data_structure() == ES_OK) {

src/core/magnetostatics/dipolar_direct_sum_gpu.cpp

@@ -38,15 +38,15 @@ static void get_simulation_box(float *box, int *per) {
 
 DipolarDirectSumGpu::DipolarDirectSumGpu(double prefactor)
     : prefactor{prefactor} {
-  auto &gpu_particle_data = System::Instance().gpu;
+  auto &gpu_particle_data = System::get_system().gpu;
   gpu_particle_data.enable_property(GpuParticleData::prop::force);
   gpu_particle_data.enable_property(GpuParticleData::prop::torque);
   gpu_particle_data.enable_property(GpuParticleData::prop::pos);
   gpu_particle_data.enable_property(GpuParticleData::prop::dip);
 }
 
 void DipolarDirectSumGpu::add_long_range_forces() const {
-  auto &gpu = System::Instance().gpu;
+  auto &gpu = System::get_system().gpu;
   gpu.update();
   if (this_node != 0) {
     return;
@@ -65,7 +65,7 @@ void DipolarDirectSumGpu::add_long_range_forces() const {
 }
 
 void DipolarDirectSumGpu::long_range_energy() const {
-  auto &gpu = System::Instance().gpu;
+  auto &gpu = System::get_system().gpu;
   gpu.update();
   if (this_node != 0) {
     return;

src/core/system/System.cpp

@@ -22,6 +22,18 @@
 
 #include <utils/Vector.hpp>
 
-System *System::m_instance = nullptr;
+#include <memory>
+
+namespace System {
+
+static std::shared_ptr<System> instance;
+
+void set_system(std::shared_ptr<System> new_instance) {
+  instance = new_instance;
+}
+
+System &get_system() { return *instance; }
 
 Utils::Vector3d System::box() const { return ::box_geo.length(); }
+
+} // namespace System

src/core/system/System.hpp

@@ -26,22 +26,21 @@
 #include <utils/Vector.hpp>
 
 #include <cstddef>
+#include <memory>
+
+namespace System {
 
 class System {
 public:
-  static System &Instance() {
-    if (!m_instance)
-      m_instance = new System;
-
-    return *m_instance;
-  }
-
 #ifdef CUDA
   GpuParticleData gpu;
 #endif // ifdef CUDA
 
   Utils::Vector3d box() const;
-
-protected:
-  static System *m_instance;
 };
+
+System &get_system();
+
+void set_system(std::shared_ptr<System> new_instance);
+
+} // namespace System

src/core/unit_tests/EspressoSystem_test.cpp

@@ -45,7 +45,7 @@
 #include <memory>
 
 namespace espresso {
-static auto &system = System::Instance();
+static auto system = ::System::System{};
 }
 
 /* Decorator to run a unit test depending on GPU availability. */