Implement replicas and dipole field tracking to DipolarDirectSumGpu

doc/sphinx/magnetostatics.rst

@@ -156,10 +156,10 @@ attach it to the system. The only required parameter is the prefactor
     dds = espressomd.magnetostatics.DipolarDirectSumGpu(prefactor=1)
     system.magnetostatics.solver = dds
 
-The CPU implementation has an optional argument ``n_replicas`` which
+Both implementations have an optional argument ``n_replicas`` which
 adds periodic copies to the system along periodic directions. In that
 case, the minimum image convention is no longer used.
-Additionally, enabling the ``DIPOLE_FIELDS_TRACKING`` feature enables the CPU
+Additionally, enabling the ``DIPOLE_FIELD_TRACKING`` feature enables either
 implementation to calculate the total dipole field at the position of each
 magnetic particle in the primary simulation box. These values are stored in
 the particle handle's ``dip_fld`` property and can be accessed directly or

src/core/energy.cpp

@@ -160,10 +160,4 @@ std::optional<double> System::particle_bond_energy(int pid, int bond_id,
   }
 }
 
-#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING
-void System::calculate_long_range_fields() {
-  dipoles.calc_long_range_field(cell_structure->local_particles());
-}
-#endif
-
 } // namespace System

src/core/forces.cpp

@@ -154,6 +154,13 @@ static void force_capping(CellStructure &cell_structure, double force_cap) {
   }
 }
 
+#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING
+static void reinit_dip_fld(CellStructure const &cell_structure) {
+  cell_structure.for_each_local_particle(
+      [](Particle &p) { p.dip_fld() = {0., 0., 0.}; });
+}
+#endif
+
 void System::System::calculate_forces() {
 #ifdef ESPRESSO_CALIPER
   CALI_CXX_MARK_FUNCTION;
@@ -182,6 +189,11 @@ void System::System::calculate_forces() {
     npt_inst_pressure->p_vir = Utils::Vector3d{};
   }
 #endif
+#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING
+  // reset dipole field
+  reinit_dip_fld(*cell_structure);
+#endif
+
   // Use combined function instead of two separate calls
 
   auto const elc_kernel = coulomb.pair_force_elc_kernel();
@@ -358,6 +370,11 @@ void System::System::calculate_forces() {
   CALI_MARK_BEGIN("copy_forces_from_GPU");
 #endif
   gpu.copy_forces_to_host(particles, this_node);
+
+#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING
+  gpu.copy_dip_fld_to_host(particles, this_node);
+#endif
+
 #ifdef ESPRESSO_CALIPER
   CALI_MARK_END("copy_forces_from_GPU");
 #endif

src/core/magnetostatics/dipolar_direct_sum.cpp

@@ -370,6 +370,9 @@ void DipolarDirectSum::add_long_range_forces(
     (*p)->force() += prefactor * fi.f;
     (*p)->torque() += prefactor * fi.torque;
   }
+#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING
+  DipolarDirectSum::dipole_field_at_part(particles);
+#endif
 }
 
 /**

src/core/magnetostatics/dipolar_direct_sum_gpu.cpp

@@ -37,8 +37,12 @@ static void get_simulation_box(BoxGeometry const &box_geo, float *box,
   }
 }
 
-DipolarDirectSumGpu::DipolarDirectSumGpu(double prefactor) {
+DipolarDirectSumGpu::DipolarDirectSumGpu(double prefactor, int n_replicas) {
   set_prefactor(prefactor);
+  this->n_replicas = n_replicas;
+  if (n_replicas < 0) {
+    throw std::domain_error("Parameter 'n_replicas' must be >= 0");
+  }
 }
 
 void DipolarDirectSumGpu::on_activation() const {
@@ -47,6 +51,9 @@ void DipolarDirectSumGpu::on_activation() const {
   gpu_particle_data.enable_property(GpuParticleData::prop::torque);
   gpu_particle_data.enable_property(GpuParticleData::prop::pos);
   gpu_particle_data.enable_property(GpuParticleData::prop::dip);
+#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING
+  gpu_particle_data.enable_property(GpuParticleData::prop::dip_fld);
+#endif
 }
 
 void DipolarDirectSumGpu::add_long_range_forces() const {
@@ -64,9 +71,14 @@ void DipolarDirectSumGpu::add_long_range_forces() const {
   auto const torques_device = gpu.get_particle_torques_device();
   auto const positions_device = gpu.get_particle_positions_device();
   auto const dipoles_device = gpu.get_particle_dipoles_device();
+  float *dipole_fields_device{nullptr};
+#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING
+  dipole_fields_device = gpu.get_particle_dip_fld_device();
+#endif
   DipolarDirectSum_kernel_wrapper_force(
       static_cast<float>(prefactor), npart, positions_device, dipoles_device,
-      forces_device, torques_device, box, periodicity);
+      dipole_fields_device, forces_device, torques_device, box, periodicity,
+      n_replicas);
 }
 
 void DipolarDirectSumGpu::long_range_energy() const {
@@ -88,4 +100,4 @@ void DipolarDirectSumGpu::long_range_energy() const {
                                          periodicity, energy_device);
 }
 
-#endif
+#endif // ESPRESSO_DIPOLAR_DIRECT_SUM

src/core/magnetostatics/dipolar_direct_sum_gpu.hpp

@@ -26,7 +26,8 @@
 #include "magnetostatics/actor.hpp"
 
 struct DipolarDirectSumGpu : public Dipoles::Actor<DipolarDirectSumGpu> {
-  DipolarDirectSumGpu(double prefactor);
+  int n_replicas;
+  DipolarDirectSumGpu(double prefactor, int n_replicas);
 
   void on_activation() const;
   void on_boxl_change() const {}