add check for wall solvers when particle centers fall outside the wall

solvers/DPStokes/mobility.h

@@ -189,6 +189,33 @@ class DPStokes : public libmobility::Mobility {
   void setPositions(device_span<const real> ipositions) override {
     this->numberParticles = ipositions.size() / 3;
     device_adapter<const real> positions(ipositions, device::cuda);
+
+    // check particles stay within walls
+    if (this->wallmode == "bottom" or this->wallmode == "slit") {
+      using namespace thrust::placeholders;
+      auto index_3 = thrust::make_transform_iterator(
+          thrust::make_counting_iterator(0), _1 * 3);
+      auto pos_z_iter =
+          thrust::make_permutation_iterator(positions.begin() + 2, index_3);
+      real min_z = thrust::reduce(
+          thrust::cuda::par, pos_z_iter, pos_z_iter + numberParticles,
+          std::numeric_limits<real>::max(), thrust::minimum<real>());
+
+      real max_z = this->dppar.zmax;
+      if (this->wallmode == "slit") {
+        max_z = thrust::reduce(
+            thrust::cuda::par, pos_z_iter, pos_z_iter + numberParticles,
+            -std::numeric_limits<real>::max(), thrust::maximum<real>());
+      }
+
+      if (min_z < this->dppar.zmin || max_z > this->dppar.zmax) {
+        throw std::runtime_error(
+            "[Mobility] The center of a particle has passed through a wall. "
+            "Check your configuration and forces to ensure particle centers "
+            "stay above the wall.");
+      }
+    }
+
     dpstokes->setPositions(positions.data(), this->numberParticles);
   }
 

solvers/NBody/mobility.h

@@ -6,10 +6,7 @@
 #include "extra/interface.h"
 #include <MobilityInterface/MobilityInterface.h>
 #include <MobilityInterface/random_finite_differences.h>
-#include <cmath>
 #include <optional>
-#include <type_traits>
-#include <vector>
 
 class NBody : public libmobility::Mobility {
   using periodicity_mode = libmobility::periodicity_mode;
@@ -117,13 +114,51 @@ class NBody : public libmobility::Mobility {
   }
 
   void setPositions(device_span<const real> ipositions) override {
-    positions.assign(ipositions.begin(), ipositions.end());
-    const auto numberParticles = this->getNumberParticles();
-    int i_Nbatch = (this->Nbatch < 0) ? 1 : this->Nbatch;
-    int i_NperBatch = (this->NperBatch < 0) ? numberParticles : this->NperBatch;
-    if (i_NperBatch * i_Nbatch != numberParticles)
+    const int numberParticles = ipositions.size() / 3;
+
+    using device_vector = thrust::device_vector<
+        real, libmobility::allocator::thrust_cached_allocator<real>>;
+
+    device_vector posZ(ipositions.begin(), ipositions.end());
+
+    if (this->kernel == nbody_rpy::kernel_type::bottom_wall) {
+
+      // shifts positions so the wall is at z=0 to match kernels
+      using namespace thrust::placeholders;
+      auto index_3 = thrust::make_transform_iterator(
+          thrust::make_counting_iterator(0), _1 * 3);
+      auto pos_z_iter =
+          thrust::make_permutation_iterator(posZ.begin() + 2, index_3);
+      if (wallHeight != 0) {
+        thrust::transform(thrust::cuda::par, pos_z_iter,
+                          pos_z_iter + numberParticles, pos_z_iter,
+                          _1 - wallHeight);
+      }
+
+      // check that all particles are above the wall
+      auto min_z = thrust::reduce(
+          thrust::cuda::par, pos_z_iter, pos_z_iter + numberParticles,
+          std::numeric_limits<real>::max(), thrust::minimum<real>());
+
+      if (min_z < 0) {
+        throw std::runtime_error(
+            "[Mobility] The center of a particle has fallen below the wall. "
+            "Check your configuration and forces to ensure particle centers "
+            "stay above the wall.");
+      }
+    }
+
+    // keep the potentially shifted positions
+    positions.assign(posZ.begin(), posZ.end());
+
+    const int i_Nbatch = (this->Nbatch < 0) ? 1 : this->Nbatch;
+    const int i_NperBatch =
+        (this->NperBatch < 0) ? numberParticles : this->NperBatch;
+
+    if (i_NperBatch * i_Nbatch != numberParticles) {
       throw std::runtime_error("[Mobility] Invalid batch parameters for NBody. "
                                "If in doubt, use the defaults.");
+    }
   }
 
   uint getNumberParticles() override { return this->positions.size() / 3; }
@@ -139,21 +174,7 @@ class NBody : public libmobility::Mobility {
           "setPositions?");
     using device_vector = thrust::device_vector<
         real, libmobility::allocator::thrust_cached_allocator<real>>;
-    device_vector posZ(positions);
-    if (wallHeight != 0) { // shifts positions so the wall is at z=0 since the
-                           // kernels are programmed as such.
-      using namespace thrust::placeholders;
-      auto index_3 = thrust::make_transform_iterator(
-          thrust::make_counting_iterator(0), _1 * 3);
-      auto iposition =
-          thrust::make_permutation_iterator(positions.begin() + 2, index_3);
-      auto opositionZ =
-          thrust::make_permutation_iterator(posZ.begin() + 2, index_3);
-      thrust::transform(thrust::cuda::par, iposition,
-                        iposition + numberParticles, opositionZ,
-                        _1 - wallHeight);
-    }
-    device_span<const real> pos(posZ);
+    device_span<const real> pos(positions);
     nbody_rpy::callBatchedNBody(pos, forces, torques, linear, angular, i_Nbatch,
                                 i_NperBatch, transMobility, rotMobility,
                                 transRotMobility, hydrodynamicRadius,

tests/test_fluctuation_dissipation.py

@@ -2,7 +2,7 @@
 import numpy as np
 from scipy.linalg import pinv, sqrtm
 from scipy.stats import kstest, norm
-from numpy.linalg import eig
+from numpy.linalg import eigh
 import logging
 
 from libMobility import PSE
@@ -29,7 +29,7 @@ def fluctuation_dissipation_KS(M, fluctuation_method):
     """
     if M.shape[0] != M.shape[1] or not np.allclose(M, M.T, rtol=0, atol=5e-5):
         raise ValueError("Matrix M must be square and symmetric.")
-    Sigma, Q = eig(M)
+    Sigma, Q = eigh(M)
     ind = np.argsort(Sigma)
     Sigma = np.sort(Sigma)
     Q = Q[:, ind]

tests/test_interface.py

@@ -1,6 +1,6 @@
 import pytest
 
-# from libMobility import *
+from libMobility import NBody, DPStokes
 import numpy as np
 from utils import (
     get_sane_params,
@@ -416,3 +416,31 @@ def test_prefactor(Solver, periodicity, includeAngular):
     assert np.allclose(divm_pf, prefactor * divm, atol=1e-4)
     if includeAngular:
         assert np.allclose(divmt_pf, prefactor * divmt, atol=1e-4)
+
+
+def test_wall_position_errors():
+    a = 1.0
+
+    solver = NBody("open", "open", "single_wall")
+    solver.setParameters(wallHeight=1.0)
+    solver.initialize(viscosity=1.0, hydrodynamicRadius=a)
+
+    pos = np.array([1.0, 2.0, 0.5])
+    with pytest.raises(RuntimeError, match="fallen below the wall"):
+        solver.setPositions(pos)
+
+    solver_one_wall = DPStokes("periodic", "periodic", "single_wall")
+    solver_one_wall.setParameters(Lx=10.0, Ly=10.0, zmin=-5.0, zmax=5.0)
+    solver_one_wall.initialize(viscosity=1.0, hydrodynamicRadius=a)
+    pos = np.array([1.0, 2.0, -5.5])
+    with pytest.raises(RuntimeError, match="passed through a wall"):
+        solver_one_wall.setPositions(pos)
+    pos = np.array([1.0, 2.0, 5.5])
+    solver_one_wall.setPositions(pos)
+
+    solver_two_wall = DPStokes("periodic", "periodic", "two_walls")
+    solver_two_wall.setParameters(Lx=10.0, Ly=10.0, zmin=-5.0, zmax=5.0)
+    solver_two_wall.initialize(viscosity=1.0, hydrodynamicRadius=a)
+    for pos in np.array([[1.0, 2.0, 5.5], [1.0, 2.0, -5.5]]):
+        with pytest.raises(RuntimeError, match="passed through a wall"):
+            solver_two_wall.setPositions(pos)