From 85e72e8f12e27a6cd8303b81f18eb369bd19c7aa Mon Sep 17 00:00:00 2001 From: Ryker Fish Date: Fri, 14 Aug 2026 17:01:04 -0600 Subject: [PATCH 1/2] add check for wall solvers when particle centers fall outside the wall --- solvers/DPStokes/mobility.h | 27 +++++++++++++++ solvers/NBody/mobility.h | 67 ++++++++++++++++++++++++------------- tests/test_interface.py | 30 ++++++++++++++++- 3 files changed, 100 insertions(+), 24 deletions(-) diff --git a/solvers/DPStokes/mobility.h b/solvers/DPStokes/mobility.h index b06b0cc0..a605649a 100644 --- a/solvers/DPStokes/mobility.h +++ b/solvers/DPStokes/mobility.h @@ -189,6 +189,33 @@ class DPStokes : public libmobility::Mobility { void setPositions(device_span ipositions) override { this->numberParticles = ipositions.size() / 3; device_adapter 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::max(), thrust::minimum()); + + 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::max(), thrust::maximum()); + } + + 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); } diff --git a/solvers/NBody/mobility.h b/solvers/NBody/mobility.h index 13cdeb2e..3cd4a42a 100644 --- a/solvers/NBody/mobility.h +++ b/solvers/NBody/mobility.h @@ -6,10 +6,7 @@ #include "extra/interface.h" #include #include -#include #include -#include -#include class NBody : public libmobility::Mobility { using periodicity_mode = libmobility::periodicity_mode; @@ -117,13 +114,51 @@ class NBody : public libmobility::Mobility { } void setPositions(device_span 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>; + + 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::max(), thrust::minimum()); + + 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>; - 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 pos(posZ); + device_span pos(positions); nbody_rpy::callBatchedNBody(pos, forces, torques, linear, angular, i_Nbatch, i_NperBatch, transMobility, rotMobility, transRotMobility, hydrodynamicRadius, diff --git a/tests/test_interface.py b/tests/test_interface.py index 3003477f..78b9bb07 100644 --- a/tests/test_interface.py +++ b/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) From 4c7eed8603e5f9bc90b1da3b02373e30d667c95b Mon Sep 17 00:00:00 2001 From: Ryker Fish Date: Mon, 17 Aug 2026 22:33:49 -0600 Subject: [PATCH 2/2] fix failing fluctuation test- eig was silently converting to a complex type causing a later function to fail. eigh explicitly ignores the imaginary parts of eigenvalues --- tests/test_fluctuation_dissipation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_fluctuation_dissipation.py b/tests/test_fluctuation_dissipation.py index 1454caf2..74c8ee5a 100644 --- a/tests/test_fluctuation_dissipation.py +++ b/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]