From 6c1f6b242b183b97f5d26288a0778ea1afff627d Mon Sep 17 00:00:00 2001 From: Ryker Fish Date: Fri, 6 Mar 2026 16:05:28 -0700 Subject: [PATCH 1/8] shrunk code for the dpstokes tests --- tests/test_dpstokes.py | 35 ++++++++++++++--------------------- 1 file changed, 14 insertions(+), 21 deletions(-) diff --git a/tests/test_dpstokes.py b/tests/test_dpstokes.py index 9961440a..00ae1a15 100644 --- a/tests/test_dpstokes.py +++ b/tests/test_dpstokes.py @@ -2,7 +2,7 @@ import libMobility as lm -def compute_with_dpstokes(pos, lx, ly, forces=None, torques=None): +def compute_with_dpstokes(pos, lx, ly, forces=None, torques=None) -> np.ndarray: params = {"viscosity": 1 / (6 * np.pi), "hydrodynamicRadius": 1.73} dpstokes = lm.DPStokes("periodic", "periodic", "single_wall") torques_on = True if torques is not None else False @@ -11,9 +11,9 @@ def compute_with_dpstokes(pos, lx, ly, forces=None, torques=None): dpstokes.setPositions(pos) mf_dp, mt_dp = dpstokes.Mdot(forces=forces, torques=torques) if torques_on: - return mt_dp + return np.array(mt_dp) else: - return mf_dp + return np.array(mf_dp) def test_non_square_box(): @@ -34,28 +34,21 @@ def test_non_square_box(): # This must be identical than doubling the box size in x and repeating the particles pos = np.vstack((pos_or, pos_or + np.array([0.5 * L_long, 0.0, 0.0]))) forces = np.vstack((forces_or, forces_or)) - mf_dp_lx = compute_with_dpstokes(pos, forces=forces, lx=L_long, ly=L_short)[ - : len(pos_or), : - ] - mt_dp_lx = compute_with_dpstokes(pos, torques=forces, lx=L_long, ly=L_short)[ - : len(pos_or), : - ] + mf_dp_lx = compute_with_dpstokes(pos, forces=forces, lx=L_long, ly=L_short) + mt_dp_lx = compute_with_dpstokes(pos, torques=forces, lx=L_long, ly=L_short) # And the same for doubling the box size in y pos = np.vstack((pos_or, pos_or + np.array([0.0, 0.5 * L_long, 0.0]))) forces = np.vstack((forces_or, forces_or)) - mf_dp_ly = compute_with_dpstokes(pos, forces=forces, lx=L_short, ly=L_long)[ - : len(pos_or), : - ] - mt_dp_ly = compute_with_dpstokes(pos, torques=forces, lx=L_short, ly=L_long)[ - : len(pos_or), : - ] - - assert np.allclose(mf_dp_lx, mf_dp_cube, rtol=1e-3, atol=1e-2) - assert np.allclose(mf_dp_ly, mf_dp_cube, rtol=1e-3, atol=1e-2) - - assert np.allclose(mt_dp_lx, mt_dp_cube, rtol=1e-3, atol=1e-2) - assert np.allclose(mt_dp_ly, mt_dp_cube, rtol=1e-3, atol=1e-2) + mf_dp_ly = compute_with_dpstokes(pos, forces=forces, lx=L_short, ly=L_long) + mt_dp_ly = compute_with_dpstokes(pos, torques=forces, lx=L_short, ly=L_long) + + slice = np.s_[0 : len(pos_or), :] + assert np.allclose(mf_dp_lx[slice], mf_dp_cube, rtol=1e-3, atol=1e-2) + assert np.allclose(mf_dp_ly[slice], mf_dp_cube, rtol=1e-3, atol=1e-2) + + assert np.allclose(mt_dp_lx[slice], mt_dp_cube, rtol=1e-3, atol=1e-2) + assert np.allclose(mt_dp_ly[slice], mt_dp_cube, rtol=1e-3, atol=1e-2) def test_isotropy(): From b3acd4069687bf0a11299b3259eb3d584553638f Mon Sep 17 00:00:00 2001 From: Ryker Fish Date: Fri, 6 Mar 2026 22:29:11 -0700 Subject: [PATCH 2/8] feat: add check for net zero forces for dpstokes --- solvers/DPStokes/extra/uammd_interface.h | 2 + solvers/DPStokes/mobility.h | 48 ++++++++++++++++++++++++ 2 files changed, 50 insertions(+) diff --git a/solvers/DPStokes/extra/uammd_interface.h b/solvers/DPStokes/extra/uammd_interface.h index 11393208..39a52633 100644 --- a/solvers/DPStokes/extra/uammd_interface.h +++ b/solvers/DPStokes/extra/uammd_interface.h @@ -3,6 +3,7 @@ */ #include #include +#include namespace uammd_dpstokes { // This is in order not to use any UAMMD related includes here. // Instead of using uammd::real I have to re define real here. @@ -41,6 +42,7 @@ struct PyParameters { // Can be either none, bottom, slit or periodic std::string mode; bool allowChangingBoxSize = false; + bool allowUnsafeForces = false; }; class DPStokesUAMMD; diff --git a/solvers/DPStokes/mobility.h b/solvers/DPStokes/mobility.h index 876179eb..52b4f4a8 100644 --- a/solvers/DPStokes/mobility.h +++ b/solvers/DPStokes/mobility.h @@ -14,6 +14,15 @@ in linear time. A. Hashemi et al. J. Chem. Phys. 158, 154101 (2023) #include #include +namespace { +struct vec3_sum { + __device__ uammd_dpstokes::real3 operator()(const uammd_dpstokes::real3 &a, + const uammd_dpstokes::real3 &b) { + return {a.x + b.x, a.y + b.y, a.z + b.z}; + } +}; +} // namespace + class DPStokes : public libmobility::Mobility { using periodicity_mode = libmobility::periodicity_mode; using Configuration = libmobility::Configuration; @@ -54,6 +63,7 @@ class DPStokes : public libmobility::Mobility { } void initialize(Parameters ipar) override { + this->par.includeAngular = ipar.includeAngular; this->dppar.viscosity = ipar.viscosity; this->lanczosTolerance = ipar.tolerance; this->dppar.mode = this->wallmode; @@ -183,6 +193,44 @@ class DPStokes : public libmobility::Mobility { device_adapter linear(ilinear, device::cuda); device_adapter angular(iangular, device::cuda); + if (this->wallmode == "nowall" && !this->dppar.allowUnsafeForces) { + real errTol = 1e-3; + bool errFlag = false; + if (iforces.size() > 0) { + const uammd_dpstokes::real3 *fptr = + reinterpret_cast(forces.data()); + uammd_dpstokes::real3 mean_f = thrust::reduce( + thrust::cuda::par, fptr, fptr + this->numberParticles, + uammd_dpstokes::real3{0, 0, 0}, vec3_sum()); + mean_f /= this->numberParticles; + if (fabs(mean_f.x) > errTol || fabs(mean_f.y) > errTol || + fabs(mean_f.z) > errTol) { + errFlag = true; + } + } + + if (this->par.includeAngular && itorques.size() > 0) { + const uammd_dpstokes::real3 *tptr = + reinterpret_cast(torques.data()); + uammd_dpstokes::real3 mean_t = thrust::reduce( + thrust::cuda::par, tptr, tptr + this->numberParticles, + uammd_dpstokes::real3{0, 0, 0}, vec3_sum()); + mean_t /= (3 * this->numberParticles); + if (fabs(mean_t.x) > errTol || fabs(mean_t.y) > errTol || + fabs(mean_t.z) > errTol) { + errFlag = true; + } + } + if (errFlag) { + throw std::runtime_error( + "[DPStokes] If using an open boundary in z, the mean force and " + "torque in each dimension must be zero or the solver will produce " + "non-physical results. If you are sure your forces and torques " + "satisfy this condition, you can set allowUnsafeForces=True to " + "bypass this check."); + } + } + dpstokes->Mdot(forces.data(), torques.data(), linear.data(), angular.data(), this->getNumberParticles(), this->getIncludeAngular()); } From 8d2b819c8764954b1bbae1c06b7e3d14880f7b11 Mon Sep 17 00:00:00 2001 From: Ryker Fish Date: Fri, 6 Mar 2026 22:29:38 -0700 Subject: [PATCH 3/8] add test for dpstokes open error checking --- tests/test_dpstokes.py | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/tests/test_dpstokes.py b/tests/test_dpstokes.py index 00ae1a15..0bf1f12c 100644 --- a/tests/test_dpstokes.py +++ b/tests/test_dpstokes.py @@ -1,5 +1,6 @@ import numpy as np import libMobility as lm +import pytest def compute_with_dpstokes(pos, lx, ly, forces=None, torques=None) -> np.ndarray: @@ -119,3 +120,21 @@ def test_non_square_matching_rpy(): results_rpy[i + 3] += mt_j assert np.allclose(results_dpstokes, results_rpy, rtol=1e-3, atol=1e-2) + + +def test_dpstokes_open_errors(): + solver = lm.DPStokes("periodic", "periodic", "open") + solver.setParameters(Lx=10.0, Ly=10.0, zmin=0.0, zmax=10.0) + solver.initialize(hydrodynamicRadius=1.0, viscosity=1.0, includeAngular=True) + + pos = np.random.uniform(0.0, 10.0, (10, 3)) + solver.setPositions(pos) + + forces = np.random.uniform(-1.0, 1.0, (10, 3)) + with pytest.raises(RuntimeError, match="DPStokes"): + solver.Mdot(forces=forces) + with pytest.raises(RuntimeError, match="DPStokes"): + solver.Mdot(torques=forces) + + forces -= np.mean(forces, axis=0) + solver.Mdot(forces=forces, torques=forces) From a520e185a34714a71bb5078ecc6e5ef415d3c65a Mon Sep 17 00:00:00 2001 From: Ryker Fish Date: Sat, 7 Mar 2026 12:34:40 -0700 Subject: [PATCH 4/8] resolve PR comments --- solvers/DPStokes/mobility.h | 19 ++++++------------- 1 file changed, 6 insertions(+), 13 deletions(-) diff --git a/solvers/DPStokes/mobility.h b/solvers/DPStokes/mobility.h index 52b4f4a8..7a13452e 100644 --- a/solvers/DPStokes/mobility.h +++ b/solvers/DPStokes/mobility.h @@ -14,15 +14,6 @@ in linear time. A. Hashemi et al. J. Chem. Phys. 158, 154101 (2023) #include #include -namespace { -struct vec3_sum { - __device__ uammd_dpstokes::real3 operator()(const uammd_dpstokes::real3 &a, - const uammd_dpstokes::real3 &b) { - return {a.x + b.x, a.y + b.y, a.z + b.z}; - } -}; -} // namespace - class DPStokes : public libmobility::Mobility { using periodicity_mode = libmobility::periodicity_mode; using Configuration = libmobility::Configuration; @@ -194,14 +185,15 @@ class DPStokes : public libmobility::Mobility { device_adapter angular(iangular, device::cuda); if (this->wallmode == "nowall" && !this->dppar.allowUnsafeForces) { - real errTol = 1e-3; + constexpr real errTol = 1e-3; bool errFlag = false; if (iforces.size() > 0) { const uammd_dpstokes::real3 *fptr = reinterpret_cast(forces.data()); uammd_dpstokes::real3 mean_f = thrust::reduce( thrust::cuda::par, fptr, fptr + this->numberParticles, - uammd_dpstokes::real3{0, 0, 0}, vec3_sum()); + uammd_dpstokes::real3{0, 0, 0}, + thrust::plus()); mean_f /= this->numberParticles; if (fabs(mean_f.x) > errTol || fabs(mean_f.y) > errTol || fabs(mean_f.z) > errTol) { @@ -214,8 +206,9 @@ class DPStokes : public libmobility::Mobility { reinterpret_cast(torques.data()); uammd_dpstokes::real3 mean_t = thrust::reduce( thrust::cuda::par, tptr, tptr + this->numberParticles, - uammd_dpstokes::real3{0, 0, 0}, vec3_sum()); - mean_t /= (3 * this->numberParticles); + uammd_dpstokes::real3{0, 0, 0}, + thrust::plus()); + mean_t /= this->numberParticles; if (fabs(mean_t.x) > errTol || fabs(mean_t.y) > errTol || fabs(mean_t.z) > errTol) { errFlag = true; From a7de6ab3324d63d67e42d37e7d272b4b53a19ae3 Mon Sep 17 00:00:00 2001 From: Ryker Fish Date: Tue, 31 Mar 2026 15:07:18 -0600 Subject: [PATCH 5/8] switch to computing sum(f) / sum(abs(f)) when calculating net force on the domain --- solvers/DPStokes/extra/uammd_interface.h | 2 + solvers/DPStokes/mobility.h | 49 ++++++++++++++++-------- 2 files changed, 34 insertions(+), 17 deletions(-) diff --git a/solvers/DPStokes/extra/uammd_interface.h b/solvers/DPStokes/extra/uammd_interface.h index 39a52633..69c4eb3b 100644 --- a/solvers/DPStokes/extra/uammd_interface.h +++ b/solvers/DPStokes/extra/uammd_interface.h @@ -10,9 +10,11 @@ namespace uammd_dpstokes { #ifndef DOUBLE_PRECISION using real = float; using real3 = float3; +using real4 = float4; #else using real = double; using real3 = double3; +using real4 = double4; #endif // This function returns either 'single' or 'double' according to the UAMMD's diff --git a/solvers/DPStokes/mobility.h b/solvers/DPStokes/mobility.h index 7a13452e..b06b0cc0 100644 --- a/solvers/DPStokes/mobility.h +++ b/solvers/DPStokes/mobility.h @@ -12,8 +12,29 @@ in linear time. A. Hashemi et al. J. Chem. Phys. 158, 154101 (2023) #include #include #include +#include #include +namespace { +struct vec3_to_vec4 { + __device__ uammd_dpstokes::real4 operator()(const uammd_dpstokes::real3 &a) { + return {a.x, a.y, a.z, abs(a.x) + abs(a.y) + abs(a.z)}; + } +}; + +uammd_dpstokes::real3 calculate_relative_vec3(const uammd_dpstokes::real3 *vptr, + uint numberParticles) { + uammd_dpstokes::real4 sum_v = thrust::transform_reduce( + thrust::cuda::par, vptr, vptr + numberParticles, vec3_to_vec4(), + uammd_dpstokes::real4{0, 0, 0, 0}, + cuda::std::plus()); + uammd_dpstokes::real3 rel_v = + uammd_dpstokes::real3{sum_v.x, sum_v.y, sum_v.z} / sum_v.w; + + return rel_v; +} +} // namespace + class DPStokes : public libmobility::Mobility { using periodicity_mode = libmobility::periodicity_mode; using Configuration = libmobility::Configuration; @@ -185,18 +206,15 @@ class DPStokes : public libmobility::Mobility { device_adapter angular(iangular, device::cuda); if (this->wallmode == "nowall" && !this->dppar.allowUnsafeForces) { - constexpr real errTol = 1e-3; + constexpr real errTol = 1e-6; bool errFlag = false; if (iforces.size() > 0) { const uammd_dpstokes::real3 *fptr = reinterpret_cast(forces.data()); - uammd_dpstokes::real3 mean_f = thrust::reduce( - thrust::cuda::par, fptr, fptr + this->numberParticles, - uammd_dpstokes::real3{0, 0, 0}, - thrust::plus()); - mean_f /= this->numberParticles; - if (fabs(mean_f.x) > errTol || fabs(mean_f.y) > errTol || - fabs(mean_f.z) > errTol) { + uammd_dpstokes::real3 rel_f = + calculate_relative_vec3(fptr, this->numberParticles); + if (fabs(rel_f.x) > errTol || fabs(rel_f.y) > errTol || + fabs(rel_f.z) > errTol) { errFlag = true; } } @@ -204,13 +222,10 @@ class DPStokes : public libmobility::Mobility { if (this->par.includeAngular && itorques.size() > 0) { const uammd_dpstokes::real3 *tptr = reinterpret_cast(torques.data()); - uammd_dpstokes::real3 mean_t = thrust::reduce( - thrust::cuda::par, tptr, tptr + this->numberParticles, - uammd_dpstokes::real3{0, 0, 0}, - thrust::plus()); - mean_t /= this->numberParticles; - if (fabs(mean_t.x) > errTol || fabs(mean_t.y) > errTol || - fabs(mean_t.z) > errTol) { + uammd_dpstokes::real3 rel_t = + calculate_relative_vec3(tptr, this->numberParticles); + if (fabs(rel_t.x) > errTol || fabs(rel_t.y) > errTol || + fabs(rel_t.z) > errTol) { errFlag = true; } } @@ -269,12 +284,12 @@ class DPStokes : public libmobility::Mobility { this->setPositions(original_pos); thrust::transform(thrust::cuda::par, thermal_drift_m.begin(), thermal_drift_m.end(), linear.begin(), linear.begin(), - thrust::plus()); + cuda::std::plus()); if (this->getIncludeAngular()) { device_adapter angular(iangular, device::cuda); thrust::transform(thrust::cuda::par, thermal_drift_d.begin(), thermal_drift_d.end(), angular.begin(), angular.begin(), - thrust::plus()); + cuda::std::plus()); } } From 1c15b60af8ee8f45c88ba5df45abdca22ad7b88d Mon Sep 17 00:00:00 2001 From: Ryker Fish Date: Tue, 31 Mar 2026 15:14:23 -0600 Subject: [PATCH 6/8] restrict cuda version --- environment.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/environment.yml b/environment.yml index 8f47f086..b79fbfe8 100644 --- a/environment.yml +++ b/environment.yml @@ -5,7 +5,7 @@ channels: dependencies: - python <3.13 - cmake >=3.24 - - cuda-version + - cuda-version<=13.1 - gxx >9,<15 - cuda-libraries-dev - cuda-nvcc From 6b6c0795785fe2a730d499d477c0da3aab49641c Mon Sep 17 00:00:00 2001 From: Ryker Fish Date: Tue, 31 Mar 2026 16:46:23 -0600 Subject: [PATCH 7/8] feat: add ability to bypass force/torque check for DPStokes open --- solvers/DPStokes/python_wrapper.cu | 10 +++++++--- tests/test_dpstokes.py | 8 ++++++++ 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/solvers/DPStokes/python_wrapper.cu b/solvers/DPStokes/python_wrapper.cu index 402e118c..c1722462 100644 --- a/solvers/DPStokes/python_wrapper.cu +++ b/solvers/DPStokes/python_wrapper.cu @@ -23,7 +23,9 @@ zmax : float allowChangingBoxSize : bool Whether the periodic extents Lx & Ly can be modified during parameter selection. Default: false. delta : float - The finite difference step size for random finite differences. Specified in units of hydrodynamicRadius. Default is 1e-3. + The finite difference step size for random finite differences. Specified in units of hydrodynamicRadius. Default: 1e-3. +allowUnsafeForces : bool + When the Z periodicity is `open` the net force/torque on the domain must be zero. Numerically, this is checked to some tolerance. Enabling this flag will by bypass the check but can lead to unphysical results if used incorrectly. Default: false. )pbdoc"; static const char *docstring = R"pbdoc( @@ -41,7 +43,7 @@ MOBILITY_PYTHONIFY_WITH_EXTRA_CODE( solver.def( "setParameters", [](DPStokes &self, real Lx, real Ly, real zmin, real zmax, - bool allowChangingBoxSize, real delta) { + bool allowChangingBoxSize, real delta, bool allowUnsafeForces) { DPStokesParameters params; params.Lx = Lx; params.Ly = Ly; @@ -49,8 +51,10 @@ MOBILITY_PYTHONIFY_WITH_EXTRA_CODE( params.zmax = zmax; params.allowChangingBoxSize = allowChangingBoxSize; params.delta = delta; + params.allowUnsafeForces = allowUnsafeForces; self.setParametersDPStokes(params); }, "Lx"_a, "Ly"_a, "zmin"_a, "zmax"_a, "allowChangingBoxSize"_a = false, - "delta"_a = 1e-3, setparameters_docstring); + "delta"_a = 1e-3, "allowUnsafeForces"_a = false, + setparameters_docstring); , docstring); diff --git a/tests/test_dpstokes.py b/tests/test_dpstokes.py index 0bf1f12c..f507b5c2 100644 --- a/tests/test_dpstokes.py +++ b/tests/test_dpstokes.py @@ -138,3 +138,11 @@ def test_dpstokes_open_errors(): forces -= np.mean(forces, axis=0) solver.Mdot(forces=forces, torques=forces) + + # check that we can bypass the error if we allow unsafe forces + forces = np.random.uniform(-1.0, 1.0, (10, 3)) + solver = lm.DPStokes("periodic", "periodic", "open") + solver.setParameters(Lx=10.0, Ly=10.0, zmin=0.0, zmax=10.0, allowUnsafeForces=True) + solver.initialize(hydrodynamicRadius=1.0, viscosity=1.0, includeAngular=True) + solver.setPositions(pos) + solver.Mdot(forces=forces, torques=forces) From 0d7cbb67f07869781708cf30891af0dbcdfa883c Mon Sep 17 00:00:00 2001 From: Ryker Fish Date: Tue, 31 Mar 2026 17:10:00 -0600 Subject: [PATCH 8/8] modify test setup to work around the new error --- tests/utils.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/tests/utils.py b/tests/utils.py index 844fad7b..510d2553 100644 --- a/tests/utils.py +++ b/tests/utils.py @@ -108,6 +108,9 @@ def get_sane_params(solverName, geom=None): params = sane_parameters["NBody_wall"].copy() else: params = sane_parameters[solverName].copy() + + if solverName == "DPStokes" and geom == "open": + params["allowUnsafeForces"] = True return params