From 2194b0a51aa76905582007d90dd79eb1c8602d1b Mon Sep 17 00:00:00 2001 From: Ryker Date: Mon, 2 Dec 2024 07:54:08 -0700 Subject: [PATCH] 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 bd4e32f4095bcc844e50f2bfe99e1bac56e6d5bb. 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 | 1 + solvers/DPStokes/extra/uammd_wrapper.cu | 1 + solvers/DPStokes/mobility.h | 1 + solvers/NBody/extra/hydrodynamicKernels.cuh | 4 + solvers/NBody/mobility.h | 1 + solvers/PSE/mobility.h | 3 + tests/test_fluctuation_dissipation.py | 12 +- tests/test_initialization.py | 11 +- tests/test_interface.py | 147 ++++++++------------ tests/test_mobility_matrix.py | 23 +-- tests/test_precision.py | 14 +- tests/test_wall_mobility.py | 27 ++-- tests/utils.py | 29 ++++ 13 files changed, 124 insertions(+), 150 deletions(-) diff --git a/solvers/DPStokes/extra/uammd_interface.h b/solvers/DPStokes/extra/uammd_interface.h index c0a7b1d..dcd5fc2 100644 --- a/solvers/DPStokes/extra/uammd_interface.h +++ b/solvers/DPStokes/extra/uammd_interface.h @@ -43,6 +43,7 @@ namespace uammd_dpstokes{ class DPStokesGlue{ std::shared_ptr dpstokes; public: + int numberParticles; //Initialize the modules with a certain set of parameters //Reinitializes if the module was already initialized diff --git a/solvers/DPStokes/extra/uammd_wrapper.cu b/solvers/DPStokes/extra/uammd_wrapper.cu index f0fcbbd..b129bb1 100644 --- a/solvers/DPStokes/extra/uammd_wrapper.cu +++ b/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 diff --git a/solvers/DPStokes/mobility.h b/solvers/DPStokes/mobility.h index e065012..14e1dc4 100644 --- a/solvers/DPStokes/mobility.h +++ b/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); } diff --git a/solvers/NBody/extra/hydrodynamicKernels.cuh b/solvers/NBody/extra/hydrodynamicKernels.cuh index dc8d342..52452a8 100644 --- a/solvers/NBody/extra/hydrodynamicKernels.cuh +++ b/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]. diff --git a/solvers/NBody/mobility.h b/solvers/NBody/mobility.h index 1355373..6601e16 100644 --- a/solvers/NBody/mobility.h +++ b/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. diff --git a/solvers/PSE/mobility.h b/solvers/PSE/mobility.h index c80e290..ade1bad 100644 --- a/solvers/PSE/mobility.h +++ b/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); } diff --git a/tests/test_fluctuation_dissipation.py b/tests/test_fluctuation_dissipation.py index 032b314..7e68d53 100644 --- a/tests/test_fluctuation_dissipation.py +++ b/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]) diff --git a/tests/test_initialization.py b/tests/test_initialization.py index c732603..a88ce42 100644 --- a/tests/test_initialization.py +++ b/tests/test_initialization.py @@ -1,7 +1,7 @@ import pytest from libMobility import * -from utils import sane_parameters +from utils import sane_parameters, solver_configs_all # List of classes inside libMobility @@ -9,14 +9,7 @@ @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) diff --git a/tests/test_interface.py b/tests/test_interface.py index 852488a..02fd12f 100644 --- a/tests/test_interface.py +++ b/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__] diff --git a/tests/test_mobility_matrix.py b/tests/test_mobility_matrix.py index 357d950..442b56b 100644 --- a/tests/test_mobility_matrix.py +++ b/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]) diff --git a/tests/test_precision.py b/tests/test_precision.py index 64ef791..05715d8 100644 --- a/tests/test_precision.py +++ b/tests/test_precision.py @@ -2,17 +2,10 @@ import numpy as np from libMobility import * - -sane_parameters = { - "PSE": {"psi": 1.0, "Lx": 32, "Ly": 32, "Lz": 32, "shearStrain": 0.0}, - "NBody": {"algorithm": "advise"}, - "SelfMobility": {"parameter": 5.0}, -} - +from utils import solver_configs_all, sane_parameters @pytest.mark.parametrize( - ("Solver", "periodicity"), - [(SelfMobility, ("open", "open", "open"))], + ("Solver", "periodicity"), solver_configs_all ) def test_precision(Solver, periodicity): # tests that Mdot gives a non-zero result when using the precision that libMobility was compiled in @@ -46,8 +39,7 @@ def test_precision(Solver, periodicity): @pytest.mark.parametrize( - ("Solver", "periodicity"), - [(SelfMobility, ("open", "open", "open"))], + ("Solver", "periodicity"), solver_configs_all ) def test_incorrect_precision(Solver, periodicity): # libMobility should work even if the inputs have an unexpected precision diff --git a/tests/test_wall_mobility.py b/tests/test_wall_mobility.py index 6b1a180..855f316 100644 --- a/tests/test_wall_mobility.py +++ b/tests/test_wall_mobility.py @@ -2,6 +2,7 @@ import numpy as np from libMobility import DPStokes, NBody +import libMobility from utils import compute_M # NOTE: Some of the following tests will only pass if compiled with double precision. @@ -12,6 +13,8 @@ "NBody": {"algorithm": "advise"}, } +precision = np.float32 if NBody.precision == "float" else np.float64 + @pytest.mark.parametrize( ("Solver", "periodicity", "tol", "ref_file"), [ @@ -21,6 +24,9 @@ ], ) def test_self_mobility_linear(Solver, periodicity, tol, ref_file): + if precision == np.float32 and Solver.__name__ == "DPStokes": + pytest.skip("The test is only valid for double precision due to how reference data was generated.") + xymax = 76.8 params = wall_params[Solver.__name__] needsTorque = False @@ -33,8 +39,6 @@ def test_self_mobility_linear(Solver, periodicity, tol, ref_file): hydrodynamicRadius = 1.0 eta = 1/4/np.sqrt(np.pi) - precision = np.float32 if Solver.precision == "float" else np.float64 - solver = Solver(*periodicity) solver.setParameters(**params) numberParticles = 1 @@ -75,6 +79,9 @@ def test_self_mobility_linear(Solver, periodicity, tol, ref_file): ], ) def test_pair_mobility_linear(Solver, periodicity, ref_file, tol): + if precision == np.float32 and Solver.__name__ == "DPStokes": + pytest.skip("The test is only valid for double precision due to how reference data was generated.") + xymax = 76.8 params = wall_params[Solver.__name__] needsTorque = False @@ -88,7 +95,6 @@ def test_pair_mobility_linear(Solver, periodicity, ref_file, tol): radH = 1.0 # hydrodynamic radius eta = 1/4/np.sqrt(np.pi) - precision = np.float32 if Solver.precision == "float" else np.float64 tol = 100*np.finfo(precision).eps solver = Solver(*periodicity) @@ -132,6 +138,9 @@ def test_pair_mobility_linear(Solver, periodicity, ref_file, tol): ], ) def test_self_mobility_angular(Solver, periodicity, ref_file): + if precision == np.float32 and Solver.__name__ == "DPStokes": + pytest.skip("The test is only valid for double precision due to how reference data was generated.") + zmax = 19.2 xymax = 76.8 params = wall_params[Solver.__name__] @@ -147,8 +156,6 @@ def test_self_mobility_angular(Solver, periodicity, ref_file): refM = ref['M'] refHeights = ref['heights'].flatten() - precision = np.float32 if Solver.precision == "float" else np.float64 - solver = Solver(*periodicity) solver.setParameters(**params) numberParticles = 1 @@ -191,6 +198,9 @@ def test_self_mobility_angular(Solver, periodicity, ref_file): ) @pytest.mark.parametrize("offset", ["x", "y"]) def test_pair_mobility_angular(Solver, periodicity, ref_file, offset): + if precision == np.float32 and Solver.__name__ == "DPStokes": + pytest.skip("The test is only valid for double precision due to how reference data was generated.") + zmax = 19.2 xymax = 76.8 params = wall_params[Solver.__name__] @@ -206,8 +216,6 @@ def test_pair_mobility_angular(Solver, periodicity, ref_file, offset): refM = ref['M'] refHeights = ref['heights'].flatten() - precision = np.float32 if Solver.precision == "float" else np.float64 - nP = 2 solver = Solver(*periodicity) solver.setParameters(**params) @@ -253,9 +261,4 @@ def test_pair_mobility_angular(Solver, periodicity, ref_file, offset): for i in range(0, nSeps): for k in range(0, nHeights): diff = abs(allM[i,k] - refM[i,k]) - temp = diff < tol - if not np.all(temp): - print(diff) - print(np.where(temp == False)) - breakpoint() assert np.all(diff < tol) \ No newline at end of file diff --git a/tests/utils.py b/tests/utils.py index b7e6134..4194896 100644 --- a/tests/utils.py +++ b/tests/utils.py @@ -1,4 +1,5 @@ import numpy as np +from libMobility import SelfMobility, PSE, NBody, DPStokes sane_parameters = { @@ -8,6 +9,34 @@ "SelfMobility": {"parameter": 5.0}, } +solver_configs_all = [ + (SelfMobility, ("open", "open", "open")), + (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")), +] + +solver_configs_torques = [ + (Solver, periodicity) + for Solver, periodicity in solver_configs_all + if not (Solver == PSE) +] + +def initialize_solver(Solver, periodicity, numberParticles, needsTorque=False): + solver = Solver(*periodicity) + solver.setParameters(**sane_parameters[Solver.__name__]) + solver.initialize( + temperature=1.0, + viscosity=1.0, + hydrodynamicRadius=1.0, + numberParticles=numberParticles, + needsTorque=needsTorque, + ) + return solver + def compute_M(solver, numberParticles, needsTorque): precision = np.float32 if solver.precision == "float" else np.float64