Torques bug fixes (#29)

* test for error when asking for angular torques without initializing correctly

* pulled out solver initialization into a function to reduce boilerplate and add test to check if positions have been set before Mdot to prevert soft failures or obscure cuda memory bug

* specify that error has to come from libmobility

* check for positions being set before running any mobility functions

* added comment (modified from a comment originally in rigidmultiblobs) about sign conventions

* Revert "check for positions being set before running any mobility functions"

This reverts commit bd4e32f.

Being done to implement a different method for preventing this error.

* added check on a per-solver basis for setPositions

* pse fix

* updated wall tests to be easier to run

* PSE fix

* simplified test parameterization

solvers/DPStokes/extra/uammd_interface.h

@@ -43,6 +43,7 @@ namespace uammd_dpstokes{
   class DPStokesGlue{
     std::shared_ptr<DPStokesUAMMD> dpstokes;
   public:
+    int numberParticles;
 
     //Initialize the modules with a certain set of parameters
     //Reinitializes if the module was already initialized

solvers/DPStokes/extra/uammd_wrapper.cu

@@ -188,6 +188,7 @@ namespace uammd_dpstokes{
   void DPStokesGlue::setPositions(const real* h_pos){
     throwIfInvalid();
     dpstokes->setPositions(h_pos);
+    this->numberParticles = dpstokes->numberParticles;
   }
 
   //Compute the dot product of the mobility matrix with the forces and/or torques acting on the previously provided positions

solvers/DPStokes/mobility.h

@@ -116,6 +116,7 @@ class DPStokes: public libmobility::Mobility{
   }
 
   void Mdot(const real* forces, const real* torques, real* linear, real* angular) override{
+    if(dpstokes->numberParticles != this->numberParticles) throw std::runtime_error("[libMobility] Wrong number of particles in positions. Did you forget to call setPositions?");
     dpstokes->Mdot(forces, torques, linear, angular);
   }
 

solvers/NBody/extra/hydrodynamicKernels.cuh

@@ -229,6 +229,8 @@ __device__ real3 RPY_WF(real3 rij, real rh){
     return correction;
   }
 
+  // note: [1] seemingly uses the left-hand rule for the cross product, so the signs are flipped on expressions
+  // in that paper that use the Levi-Civita symbol
   __device__ real3 wallCorrection_UT(real3 rij, bool self, real h, real3 vj){
       real3 correction = real3();
       if(self){ // B2^T*vj in [1]. ^T denotes transpose.
@@ -260,6 +262,8 @@ __device__ real3 RPY_WF(real3 rij, real rh){
   return correction;
 }
 
+  // note: [1] seemingly uses the left-hand rule for the cross product, so the signs are flipped on expressions
+  // in that paper that use the Levi-Civita symbol
   __device__ real3 wallCorrection_WF(real3 rij, bool self, real h, real3 vj){
         real3 correction = real3();
         if(self){ // B2*fj in [1].

solvers/NBody/mobility.h

@@ -78,6 +78,7 @@ class NBody: public libmobility::Mobility{
 
   virtual void Mdot(const real* forces, const real* torques, real* linear, real* angular) override{
     int numberParticles = positions.size()/3;
+    if(numberParticles != this->numberParticles) throw std::runtime_error("[libMobility] Wrong number of particles in positions. Did you forget to call setPositions?");
     // Donev: Why can't there be a single callBatchedNBody routine that does if(kernel == kernel_type::bottom_wall) internally?
     // Example, what if in the future we add manually periodized RPY where one repeats a unit cell a certain number of times in each direction. This is actually easy to do with 3 loops and useful, and only requires adding a parameter nUnitCellsRepeat[3] and removing the error if some direction is periodic and only spitting an error if two walls are asked for.
 

solvers/PSE/mobility.h

@@ -74,12 +74,14 @@ class PSE: public libmobility::Mobility{
 		    real* linear, real* angular) override{
     if(torques)
       throw std::runtime_error("[PSE] Torque is not implemented");
+    if(positions.size() != 3*currentNumberParticles) throw std::runtime_error("[libMobility] Wrong number of particles in positions. Did you forget to call setPositions?");
     pse->computeHydrodynamicDisplacements(positions.data(), forces, linear, 0, 0);
   }
 
   virtual void sqrtMdotW(real* linear, real *angular, real prefactor = 1) override{
     if(angular)
       throw std::runtime_error("[PSE] Torque is not implemented");
+    if(positions.size() != 3*currentNumberParticles) throw std::runtime_error("[libMobility] Wrong number of particles in positions. Did you forget to call setPositions?");
     pse->computeHydrodynamicDisplacements(positions.data(), nullptr, linear,
 					  temperature, prefactor);
   }
@@ -91,6 +93,7 @@ class PSE: public libmobility::Mobility{
 				      real prefactor = 1) override{
     if(angular)
       throw std::runtime_error("[PSE] Torque is not implemented");
+    if(positions.size() != 3*currentNumberParticles) throw std::runtime_error("[libMobility] Wrong number of particles in positions. Did you forget to call setPositions?");
     pse->computeHydrodynamicDisplacements(positions.data(), forces, linear,
 					  temperature, prefactor);
   }

tests/test_fluctuation_dissipation.py

@@ -6,7 +6,7 @@
 import logging
 
 from libMobility import SelfMobility, PSE, NBody, DPStokes
-from utils import compute_M, sane_parameters, generate_positions_in_box
+from utils import compute_M, sane_parameters, generate_positions_in_box, solver_configs_all
 
 
 def fluctuation_dissipation_KS(M, fluctuation_method, needsTorques):
@@ -58,15 +58,7 @@ def fluctuation_dissipation_KS(M, fluctuation_method, needsTorques):
 
 
 @pytest.mark.parametrize(
-    ("Solver", "periodicity"),
-    [
-        (SelfMobility, ("open", "open", "open")),
-        (PSE, ("periodic", "periodic", "periodic")),
-        (NBody, ("open", "open", "open")),
-        (DPStokes, ("periodic", "periodic", "open")),
-        (DPStokes, ("periodic", "periodic", "single_wall")),
-        (DPStokes, ("periodic", "periodic", "two_walls")),
-    ],
+    ("Solver", "periodicity"), solver_configs_all
 )
 @pytest.mark.parametrize("hydrodynamicRadius", [0.95, 1.12])
 @pytest.mark.parametrize("numberParticles", [1,2,10])

tests/test_initialization.py

@@ -1,22 +1,15 @@
 import pytest
 
 from libMobility import *
-from utils import sane_parameters
+from utils import sane_parameters, solver_configs_all
 
 # List of classes inside libMobility
 
 solver_list = [PSE, NBody, DPStokes, SelfMobility]
 
 
 @pytest.mark.parametrize(
-    ("Solver", "periodicity"),
-    [
-        (PSE, ("periodic", "periodic", "periodic")),
-        (NBody, ("open", "open", "open")),
-        (DPStokes, ("periodic", "periodic", "open")),
-        (DPStokes, ("periodic", "periodic", "single_wall")),
-        (DPStokes, ("periodic", "periodic", "two_walls")),
-    ],
+    ("Solver", "periodicity"), solver_configs_all
 )
 def test_periodic_initialization(Solver, periodicity):
     solver = Solver(*periodicity)

tests/test_interface.py

@@ -1,32 +1,16 @@
 import pytest
 from libMobility import *
 import numpy as np
-from utils import sane_parameters
+from utils import sane_parameters, initialize_solver, solver_configs_all, solver_configs_torques
 
 
 @pytest.mark.parametrize(
-    ("Solver", "periodicity"),
-    [
-        (SelfMobility, ("open", "open", "open")),
-        (PSE, ("periodic", "periodic", "periodic")),
-        (NBody, ("open", "open", "open")),
-        (DPStokes, ("periodic", "periodic", "open")),
-        (DPStokes, ("periodic", "periodic", "single_wall")),
-        (DPStokes, ("periodic", "periodic", "two_walls")),
-    ],
+    ("Solver", "periodicity"), solver_configs_all
 )
 def test_contiguous(Solver, periodicity):
-    hydrodynamicRadius = 1.0
-    solver = Solver(*periodicity)
-    parameters = sane_parameters[Solver.__name__]
-    solver.setParameters(**parameters)
+
     numberParticles = 1
-    solver.initialize(
-        temperature=1.0,
-        viscosity=1.0,
-        hydrodynamicRadius=hydrodynamicRadius,
-        numberParticles=numberParticles,
-    )
+    solver = initialize_solver(Solver, periodicity, numberParticles)
 
     # Set precision to be the same as compiled precision
     precision = np.float32 if Solver.precision == "float" else np.float64
@@ -40,28 +24,12 @@ def test_contiguous(Solver, periodicity):
 
 
 @pytest.mark.parametrize(
-    ("Solver", "periodicity"),
-    [
-        (SelfMobility, ("open", "open", "open")),
-        (PSE, ("periodic", "periodic", "periodic")),
-        (NBody, ("open", "open", "open")),
-        (DPStokes, ("periodic", "periodic", "open")),
-        (DPStokes, ("periodic", "periodic", "single_wall")),
-        (DPStokes, ("periodic", "periodic", "two_walls")),
-    ],
+    ("Solver", "periodicity"), solver_configs_all
 )
 def test_returns_mf(Solver, periodicity):
-    hydrodynamicRadius = 1.0
-    solver = Solver(*periodicity)
-    parameters = sane_parameters[Solver.__name__]
-    solver.setParameters(**parameters)
+
     numberParticles = 1
-    solver.initialize(
-        temperature=1.0,
-        viscosity=1.0,
-        hydrodynamicRadius=hydrodynamicRadius,
-        numberParticles=numberParticles,
-    )
+    solver = initialize_solver(Solver, periodicity, numberParticles)
 
     # Set precision to be the same as compiled precision
     precision = np.float32 if Solver.precision == "float" else np.float64
@@ -76,28 +44,12 @@ def test_returns_mf(Solver, periodicity):
 
 
 @pytest.mark.parametrize(
-    ("Solver", "periodicity"),
-    [
-        (SelfMobility, ("open", "open", "open")),
-        (NBody, ("open", "open", "open")),
-        (DPStokes, ("periodic", "periodic", "open")),
-        (DPStokes, ("periodic", "periodic", "single_wall")),
-        (DPStokes, ("periodic", "periodic", "two_walls")),
-    ],
+    ("Solver", "periodicity"), solver_configs_torques
 )
 def test_returns_mf_mt(Solver, periodicity):
-    hydrodynamicRadius = 1.0
-    solver = Solver(*periodicity)
-    parameters = sane_parameters[Solver.__name__]
-    solver.setParameters(**parameters)
+
     numberParticles = 1
-    solver.initialize(
-        temperature=1.0,
-        viscosity=1.0,
-        hydrodynamicRadius=hydrodynamicRadius,
-        numberParticles=numberParticles,
-        needsTorque=True,
-    )
+    solver = initialize_solver(Solver, periodicity, numberParticles, needsTorque=True)
 
     # Set precision to be the same as compiled precision
     precision = np.float32 if Solver.precision == "float" else np.float64
@@ -116,28 +68,12 @@ def test_returns_mf_mt(Solver, periodicity):
 
 
 @pytest.mark.parametrize(
-    ("Solver", "periodicity"),
-    [
-        (SelfMobility, ("open", "open", "open")),
-        (PSE, ("periodic", "periodic", "periodic")),
-        (NBody, ("open", "open", "open")),
-        (DPStokes, ("periodic", "periodic", "open")),
-        (DPStokes, ("periodic", "periodic", "single_wall")),
-        (DPStokes, ("periodic", "periodic", "two_walls")),
-    ],
+    ("Solver", "periodicity"), solver_configs_all
 )
 def test_returns_sqrtM(Solver, periodicity):
-    hydrodynamicRadius = 1.0
-    solver = Solver(*periodicity)
-    parameters = sane_parameters[Solver.__name__]
-    solver.setParameters(**parameters)
+
     numberParticles = 1
-    solver.initialize(
-        temperature=1.0,
-        viscosity=1.0,
-        hydrodynamicRadius=hydrodynamicRadius,
-        numberParticles=numberParticles,
-    )
+    solver = initialize_solver(Solver, periodicity, numberParticles)
 
     # Set precision to be the same as compiled precision
     precision = np.float32 if Solver.precision == "float" else np.float64
@@ -148,15 +84,7 @@ def test_returns_sqrtM(Solver, periodicity):
 
 
 @pytest.mark.parametrize(
-    ("Solver", "periodicity"),
-    [
-        (SelfMobility, ("open", "open", "open")),
-        (PSE, ("periodic", "periodic", "periodic")),
-        (NBody, ("open", "open", "open")),
-        (DPStokes, ("periodic", "periodic", "open")),
-        (DPStokes, ("periodic", "periodic", "single_wall")),
-        (DPStokes, ("periodic", "periodic", "two_walls")),
-    ],
+    ("Solver", "periodicity"), solver_configs_all
 )
 def test_returns_hydrodisp(Solver, periodicity):
     hydrodynamicRadius = 1.0
@@ -181,21 +109,64 @@ def test_returns_hydrodisp(Solver, periodicity):
     sqrtmw, _ = solver.hydrodynamicVelocities(forces)
     assert sqrtmw.shape == (numberParticles, 3)
 
+@pytest.mark.parametrize(
+    ("Solver", "periodicity"), solver_configs_torques
+)
+def test_no_torques_error(Solver, periodicity):
+    # Test that the solver raises an error if torques are provided but solver was not initialized with needsTorque=True
+    numberParticles = 1
+    solver = initialize_solver(Solver, periodicity, numberParticles, needsTorque=False)
+
+    # Set precision to be the same as compiled precision
+    precision = np.float32 if Solver.precision == "float" else np.float64
+
+    forces = np.random.rand(numberParticles, 3).astype(precision)
+    torques = np.random.rand(numberParticles, 3).astype(precision)
+    posistions = np.random.rand(numberParticles, 3).astype(precision)
+
+    solver.setPositions(posistions)
+
+    with pytest.raises(RuntimeError):
+        u, w = solver.Mdot(forces, torques)
 
 @pytest.mark.parametrize(
     ("Solver", "periodicity"),
     [
-        (SelfMobility, ("open", "open", "open")),
-        (PSE, ("periodic", "periodic", "periodic")),
         (NBody, ("open", "open", "open")),
+        (NBody, ("open", "open", "single_wall")),
+        (PSE, ("periodic", "periodic", "periodic")),
         (DPStokes, ("periodic", "periodic", "open")),
         (DPStokes, ("periodic", "periodic", "single_wall")),
         (DPStokes, ("periodic", "periodic", "two_walls")),
     ],
 )
+def test_no_positions_error(Solver, periodicity):
+    # Test that the solver raises an error if Mdot, sqrtMdot, and hydroDisp are called before setting positions
+    numberParticles = 1
+    solver = initialize_solver(Solver, periodicity, numberParticles)
+
+    # Set precision to be the same as compiled precision
+    precision = np.float32 if Solver.precision == "float" else np.float64
+
+    forces = np.random.rand(numberParticles, 3).astype(precision)
+
+    with pytest.raises(RuntimeError, match=r"\[libMobility\]*"):
+        u, _ = solver.Mdot(forces)
+
+    with pytest.raises(RuntimeError, match=r"\[libMobility\]*"):
+        sqrtmw, _ = solver.sqrtMdotW()
+
+    with pytest.raises(RuntimeError, match=r"\[libMobility\]*"):
+        u, _ = solver.hydrodynamicVelocities()
+
+@pytest.mark.parametrize(
+    ("Solver", "periodicity"), solver_configs_all
+)
 @pytest.mark.parametrize("needsTorque", [True, False])
 def test_hydrodisp_equivalent(Solver, periodicity, needsTorque):
     #  Check that calling Mdot is equivalent to calling hydrodynamicVelocities with temperature = 0
+    if needsTorque and Solver.__name__ == "PSE":
+        pytest.skip("PSE does not support torques")
     hydrodynamicRadius = 1.0
     solver = Solver(*periodicity)
     parameters = sane_parameters[Solver.__name__]

tests/test_mobility_matrix.py

@@ -2,20 +2,11 @@
 import pytest
 import numpy as np
 from libMobility import SelfMobility, PSE, NBody, DPStokes
-from utils import compute_M
+from utils import compute_M, solver_configs_all, solver_configs_torques
 
 
 @pytest.mark.parametrize(
-    ("Solver", "periodicity"),
-    [
-        (SelfMobility, ("open", "open", "open")),
-        (PSE, ("periodic", "periodic", "periodic")),
-        (NBody, ("open", "open", "open")),
-        (NBody, ("open", "open", "single_wall")),
-        (DPStokes, ("periodic", "periodic", "open")),
-        (DPStokes, ("periodic", "periodic", "single_wall")),
-        (DPStokes, ("periodic", "periodic", "two_walls")),
-    ],
+    ("Solver", "periodicity"), solver_configs_torques
 )
 @pytest.mark.parametrize("hydrodynamicRadius", [1.0, 0.95, 1.12])
 @pytest.mark.parametrize("numberParticles", [1, 2, 3, 10])
@@ -48,15 +39,7 @@ def test_mobility_matrix_linear(
 
 
 @pytest.mark.parametrize(
-    ("Solver", "periodicity"),
-    [
-        (SelfMobility, ("open", "open", "open")),
-        (NBody, ("open", "open", "open")),
-        (NBody, ("open", "open", "single_wall")),
-        (DPStokes, ("periodic", "periodic", "open")),
-        (DPStokes, ("periodic", "periodic", "single_wall")),
-        (DPStokes, ("periodic", "periodic", "two_walls")),
-    ],
+    ("Solver", "periodicity"), solver_configs_torques
 )
 @pytest.mark.parametrize("hydrodynamicRadius", [1.0, 0.95, 1.12])
 @pytest.mark.parametrize("numberParticles", [1, 2, 3, 10])