Added check if particles are not in ELC gap region (#4051)

Fixes Issue #4009

Description of changes:
- add a check in both positive and negative direction with error message

doc/sphinx/electrostatics.rst

@@ -255,8 +255,8 @@ Usage notes:
     z-direction not to contain particles. The size in z-direction of this slab
     is controlled by the ``gap_size`` parameter. The user has to ensure that
     no particles enter this region by means of constraints or by fixing the
-    particles' z-coordinate. When there is no empty slab of the specified size,
-    the method will silently produce wrong results.
+    particles' z-coordinate. When particles enter the slab of the specified
+    size, an error will be thrown.
 
 *ELC* is an |es| actor and is used with::
 

doc/sphinx/magnetostatics.rst

@@ -96,8 +96,8 @@ Usage notes:
     z-direction not to contain particles. The size in z-direction of this slab
     is controlled by the ``gap_size`` parameter. The user has to ensure that
     no particles enter this region by means of constraints or by fixing the
-    particles' z-coordinate. When there is no empty slab of the specified size,
-    the method will silently produce wrong results.
+    particles' z-coordinate. When particles enter the slab of the specified
+    size, an error will be thrown.
 
   * The method can be tuned using the ``accuracy`` parameter. In contrast to
     the electrostatic method, it refers to the energy. Furthermore, it is

src/core/electrostatics_magnetostatics/elc.cpp

@@ -223,6 +223,20 @@ void distribute(int size) {
   MPI_Allreduce(send_buf, gblcblk, size, MPI_DOUBLE, MPI_SUM, comm_cart);
 }
 
+/** Checks if a charged particle is in the forbidden gap region
+ */
+inline void check_gap_elc(const Particle &p) {
+  if (p.p.q != 0) {
+    if (p.r.p[2] < 0)
+      runtimeErrorMsg() << "Particle " << p.p.identity << " entered ELC gap "
+                        << "region by " << (p.r.p[2]);
+    else if (p.r.p[2] > elc_params.h) {
+      runtimeErrorMsg() << "Particle " << p.p.identity << " entered ELC gap "
+                        << "region by " << (p.r.p[2] - elc_params.h);
+    }
+  }
+}
+
 /*****************************************************************/
 /* dipole terms */
 /*****************************************************************/
@@ -247,6 +261,8 @@ static void add_dipole_force(const ParticleRange &particles) {
   gblcblk[2] = 0; // sum q_i
 
   for (auto const &p : local_particles) {
+    check_gap_elc(p);
+
     gblcblk[0] += p.p.q * (p.r.p[2] - shift);
     gblcblk[1] += p.p.q * p.r.p[2];
     gblcblk[2] += p.p.q;
@@ -311,6 +327,8 @@ static double dipole_energy(const ParticleRange &particles) {
   gblcblk[6] = 0; // sum q_i z_i           primary box
 
   for (auto &p : particles) {
+    check_gap_elc(p);
+
     gblcblk[0] += p.p.q;
     gblcblk[2] += p.p.q * (p.r.p[2] - shift);
     gblcblk[4] += p.p.q * (Utils::sqr(p.r.p[2] - shift));

src/core/electrostatics_magnetostatics/mdlc_correction.cpp

@@ -42,6 +42,20 @@
 
 DLC_struct dlc_params = {1e100, 0, 0, 0, 0};
 
+/** Checks if a magnetic particle is in the forbidden gap region
+ */
+inline void check_gap_mdlc(const Particle &p) {
+  if (p.p.dipm != 0) {
+    if (p.r.p[2] < 0)
+      runtimeErrorMsg() << "Particle " << p.p.identity << " entered MDLC gap "
+                        << "region by " << (p.r.p[2]);
+    else if (p.r.p[2] > dlc_params.h) {
+      runtimeErrorMsg() << "Particle " << p.p.identity << " entered MDLC gap "
+                        << "region by " << (p.r.p[2] - dlc_params.h);
+    }
+  }
+}
+
 /** Calculate the maximal dipole moment in the system */
 double calc_mu_max() {
   auto const local_particles = cell_structure.local_particles();
@@ -341,6 +355,8 @@ void add_mdlc_force_corrections(const ParticleRange &particles) {
 
   int ip = 0;
   for (auto &p : particles) {
+    check_gap_mdlc(p);
+
     if ((p.p.dipm) != 0.0) {
       // SDC correction term is zero for the forces
       p.f.f += dipole.prefactor * dip_DLC_f[ip];
@@ -370,6 +386,13 @@ double add_mdlc_energy_corrections(const ParticleRange &particles) {
 
   auto const volume = box_geo.volume();
 
+  // Check if particles aren't in the forbidden gap region
+  // This loop is needed, because there is no other guaranteed
+  // single pass over all particles in this function.
+  for (auto const &p : particles) {
+    check_gap_mdlc(p);
+  }
+
   //---- Compute the corrections ----------------------------------
 
   // First the DLC correction

testsuite/python/dipolar_mdlc_p3m_scafacos_p2nfft.py

@@ -62,6 +62,7 @@ def test_mdlc(self):
         s.box_l = 3 * [box_l]
         ref_E_path = abspath("data/mdlc_reference_data_energy.dat")
         ref_E = float(np.genfromtxt(ref_E_path))
+        gap_size = 2.0
 
         # Particles
         data = np.genfromtxt(
@@ -70,7 +71,7 @@ def test_mdlc(self):
         s.part[:].rotation = (1, 1, 1)
 
         p3m = magnetostatics.DipolarP3M(prefactor=1, mesh=32, accuracy=1E-4)
-        dlc = magnetostatic_extensions.DLC(maxPWerror=1E-5, gap_size=2.)
+        dlc = magnetostatic_extensions.DLC(maxPWerror=1E-5, gap_size=gap_size)
         s.actors.add(p3m)
         s.actors.add(dlc)
         s.integrator.run(0)
@@ -86,6 +87,23 @@ def test_mdlc(self):
         self.assertLessEqual(abs(err_t), tol_t, "Torque difference too large")
         self.assertLessEqual(abs(err_f), tol_f, "Force difference too large")
 
+        # Check if error is thrown when particles enter the MDLC gap
+        # positive direction
+        s.part[0].pos = [
+            s.box_l[0] / 2,
+            s.box_l[1] / 2,
+            s.box_l[2] - gap_size / 2]
+        with self.assertRaises(Exception):
+            self.system.analysis.energy()
+        with self.assertRaises(Exception):
+            self.integrator.run(2)
+        # negative direction
+        s.part[0].pos = [s.box_l[0] / 2, s.box_l[1] / 2, -gap_size / 2]
+        with self.assertRaises(Exception):
+            self.system.analysis.energy()
+        with self.assertRaises(Exception):
+            self.integrator.run(2)
+
     def test_p3m(self):
         s = self.system
         rho = 0.09

testsuite/python/elc.py

@@ -71,6 +71,20 @@ def test_finite_potential_drop(self):
         self.assertAlmostEqual(E_expected, p1.f[2] / p1.q)
         self.assertAlmostEqual(E_expected, p2.f[2] / p2.q)
 
+        # Check if error is thrown when particles enter the ELC gap
+        # positive direction
+        p1.pos = [BOX_L[0] / 2, BOX_L[1] / 2, BOX_L[2] - GAP[2] / 2]
+        with self.assertRaises(Exception):
+            self.system.analysis.energy()
+        with self.assertRaises(Exception):
+            self.integrator.run(2)
+        # negative direction
+        p1.pos = [BOX_L[0] / 2, BOX_L[1] / 2, -GAP[2] / 2]
+        with self.assertRaises(Exception):
+            self.system.analysis.energy()
+        with self.assertRaises(Exception):
+            self.integrator.run(2)
+
 
 if __name__ == "__main__":
     ut.main()