From e77c7c3585c9e4e16903ab877ccdf390bd1f9c69 Mon Sep 17 00:00:00 2001 From: Raul Date: Sat, 4 May 2024 23:51:09 +0200 Subject: [PATCH] Make the different functions return a numpy array instead of requiring it as input (#15) * Update CMakeLists.txt * Move examples to a new folder. Change the location of some installation headers * Update CI * Update READMEs * Update README * Add option for cuda archs * Update Python interface * Add a function to get the number of particles * Change the interface so that Mdot returns an array * Update docs * Port all functions to return results. Add some tests for the interface Do not call Mdot if forces is null * Update env * Relax python dependency * Add LICENSE * Update test * Fix test * Make fluctuation test more robust --- CMakeLists.txt | 5 +- LICENSE | 7 + docs/source/usage.rst | 5 +- environment.yml | 3 +- include/MobilityInterface/MobilityInterface.h | 8 +- include/MobilityInterface/pythonify.h | 241 +++++++++++------- solvers/SelfMobility/mobility.h | 3 + tests/test_fluctuation_dissipation.py | 6 +- tests/test_interface.py | 120 ++++++++- tests/test_mobility_matrix.py | 9 +- tests/test_precision.py | 12 +- tests/utils.py | 4 +- 12 files changed, 300 insertions(+), 123 deletions(-) create mode 100644 LICENSE diff --git a/CMakeLists.txt b/CMakeLists.txt index 04e8a4f..efe45a0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,7 +3,6 @@ project(libMobility) enable_language(CUDA) set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/" "$ENV{CONDA_PREFIX}/share/cmake/Modules") set(CMAKE_BUILD_TYPE Release) - add_compile_definitions(PUBLIC MAXLOGLEVEL=3) set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_STANDARD 14) @@ -11,7 +10,9 @@ set(CMAKE_CUDA_STANDARD 14) set(CMAKE_CUDA_STANDARD_REQUIRED ON) set(CMAKE_CUDA_SEPARABLE_COMPILATION OFF) # Set CUDA archs so all supported GPUs are covered -set(CMAKE_CUDA_ARCHITECTURES "all") +if(NOT DEFINED CMAKE_CUDA_ARCHITECTURES) + set(CMAKE_CUDA_ARCHITECTURES "all") +endif() #add_compile_definitions(PUBLIC DOUBLE_PRECISION) set(UAMMD_INCLUDE include/third_party/uammd/src include/third_party/uammd/src/third_party) diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..83702c8 --- /dev/null +++ b/LICENSE @@ -0,0 +1,7 @@ +Copyright 2024 StochasticHydroTools + +Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. diff --git a/docs/source/usage.rst b/docs/source/usage.rst index d100aa5..dc71d80 100644 --- a/docs/source/usage.rst +++ b/docs/source/usage.rst @@ -62,7 +62,6 @@ A libMobility solver is initialized in three steps: precision = np.float32 if Solver.precision=="float" else np.float64 pos = np.random.rand(3*numberParticles).astype(precision) force = np.ones(3*numberParticles).astype(precision) - result = np.zeros(3*numberParticles).astype(precision) # The solver will fail if it is not compatible with the provided periodicity nb = Solver(periodicityX='open',periodicityY='open',periodicityZ='open') @@ -72,11 +71,11 @@ A libMobility solver is initialized in three steps: hydrodynamicRadius = 1.0, numberParticles = numberParticles) nb.setPositions(pos) - nb.Mdot(forces = force, result = result) + result = nb.Mdot(forces = force) print(f"{numberParticles} particles located at ( X Y Z ): {pos}") print("Forces:", force) print("M*F:", result) # result = prefactor*sqrt(2*temperature*M)*dW - nb.sqrtMdotW(prefactor = 1.0, result = result) + result = nb.sqrtMdotW(prefactor = 1.0) print("sqrt(2*T*M)*N(0,1):", result) nb.clean() diff --git a/environment.yml b/environment.yml index 6aff019..c3bd25d 100644 --- a/environment.yml +++ b/environment.yml @@ -2,7 +2,7 @@ name: libmobility channels: - conda-forge dependencies: - - cmake >=3.22 + - cmake >=3.24 - cuda-version 12.* - gxx - cuda-libraries-dev @@ -11,6 +11,5 @@ dependencies: - mkl-devel - pybind11 - numpy - - python 3.11.* - pytest - scipy diff --git a/include/MobilityInterface/MobilityInterface.h b/include/MobilityInterface/MobilityInterface.h index 969ad35..29e7d7c 100644 --- a/include/MobilityInterface/MobilityInterface.h +++ b/include/MobilityInterface/MobilityInterface.h @@ -107,10 +107,16 @@ namespace libmobility{ //Equivalent to calling Mdot and then stochasticDisplacements, can be faster in some solvers virtual void hydrodynamicVelocities(const real* forces, real* result, real prefactor = 1){ - Mdot(forces, result); + if(forces){ + Mdot(forces, result); + } sqrtMdotW(result, prefactor); } + int getNumberParticles() const{ + return this->numberParticles; + } + //Clean any memory allocated by the solver virtual void clean(){ lanczos.reset(); diff --git a/include/MobilityInterface/pythonify.h b/include/MobilityInterface/pythonify.h index 59465b2..653ed06 100644 --- a/include/MobilityInterface/pythonify.h +++ b/include/MobilityInterface/pythonify.h @@ -1,29 +1,41 @@ /*Raul P. Pelaez 2021. -The MOBILITY_PYTHONIFY(className, description) macro creates a pybind11 module from a class (called className) that inherits from libmobility::Mobility. "description" is a string that will be printed when calling help(className) from python (accompanied by the default documentation of the mobility interface. +The MOBILITY_PYTHONIFY(className, description) macro creates a pybind11 module +from a class (called className) that inherits from libmobility::Mobility. +"description" is a string that will be printed when calling help(className) from +python (accompanied by the default documentation of the mobility interface. */ #include #ifndef MOBILITY_PYTHONIFY_H -#include"MobilityInterface.h" -#include -#include +#include "MobilityInterface.h" +#include +#include namespace py = pybind11; using namespace pybind11::literals; using pyarray = py::array; +using pyarray_c = + py::array_t; #define MOBILITYSTR(s) xMOBILITYSTR(s) #define xMOBILITYSTR(s) #s -inline auto string2Periodicity(std::string per){ +inline auto string2Periodicity(std::string per) { using libmobility::periodicity_mode; - if(per == "open") return periodicity_mode::open; - else if(per == "unspecified") return periodicity_mode::unspecified; - else if(per == "single_wall") return periodicity_mode::single_wall; - else if(per == "two_walls") return periodicity_mode::two_walls; - else if(per == "periodic") return periodicity_mode::periodic; - else throw std::runtime_error("[libMobility] Invalid periodicity"); + if (per == "open") + return periodicity_mode::open; + else if (per == "unspecified") + return periodicity_mode::unspecified; + else if (per == "single_wall") + return periodicity_mode::single_wall; + else if (per == "two_walls") + return periodicity_mode::two_walls; + else if (per == "periodic") + return periodicity_mode::periodic; + else + throw std::runtime_error("[libMobility] Invalid periodicity"); } -inline auto createConfiguration(std::string perx, std::string pery, std::string perz){ +inline auto createConfiguration(std::string perx, std::string pery, + std::string perz) { libmobility::Configuration conf; conf.periodicityX = string2Periodicity(perx); conf.periodicityY = string2Periodicity(pery); @@ -31,26 +43,27 @@ inline auto createConfiguration(std::string perx, std::string pery, std::string return conf; } -template -void check_dtype(pyarray &arr){ - if(not py::isinstance>(arr)){ +template void check_dtype(Array &arr) { + if (not py::isinstance>(arr)) { throw std::runtime_error("Input array must have the correct data type."); } - if(not py::isinstance>(arr)){ - throw std::runtime_error("The input array is not contiguous and cannot be used as a buffer."); + if (not py::isinstance< + py::array_t>(arr)) { + throw std::runtime_error( + "The input array is not contiguous and cannot be used as a buffer."); } - } -libmobility::real* cast_to_real(pyarray &arr){ +template libmobility::real *cast_to_real(Array &arr) { check_dtype(arr); - return static_cast(arr.mutable_data()); + return static_cast(arr.mutable_data()); } -const libmobility::real* cast_to_const_real(pyarray &arr){ +template const libmobility::real *cast_to_const_real(Array &arr) { check_dtype(arr); - return static_cast(arr.data()); + return static_cast(arr.data()); } + const char *constructor_docstring = R"pbdoc( Initialize the module with a given set of periodicity conditions. @@ -88,39 +101,90 @@ tolerance : float, optional Tolerance, used for approximate methods and also for Lanczos (default fluctuation computation). Default is 1e-4. )pbdoc"; -const char *mdot_docstring = R"pbdoc( -Computes the product of the Mobility matrix with a group of forces, :math:`\boldsymbol{\mathcal{M}}\boldsymbol{F}`. +template +auto call_sqrtMdotW(Solver &solver, libmobility::real prefactor) { + auto result = + py::array_t({solver.getNumberParticles() * 3}); + solver.sqrtMdotW(cast_to_real(result), prefactor); + return result.reshape({solver.getNumberParticles(), 3}); +} + +const char *sqrtMdotW_docstring = R"pbdoc( +Computes the stochastic contribution, :math:`\text{prefactor}\sqrt{2T\boldsymbol{\mathcal{M}}}d\boldsymbol{W}`, where :math:`\boldsymbol{\mathcal{M}}` is the mobility matrix and :math:`d\boldsymbol{W}` is a Wiener process. It is required that :py:mod:`setPositions` has been called before calling this function. -Both inputs must have precision given by the precision attribute of the module. -Both inputs must have size 3*N, where N is the number of particles. -The arrays are ordered as :code:`[f0x, f0y, f0z, f1x, f1y, f1z, ...]`. Parameters ---------- -forces : array_like, - Forces acting on the particles. -result : array_like, - Where the result will be stored. + +prefactor : float, optional + Prefactor to multiply the result by. Default is 1.0. + +Returns +------- +array_like + The resulting fluctuations. Shape is (N, 3), where N is the number of particles. + )pbdoc"; -const char *sqrtMdotW_docstring = R"pbdoc( -Computes the stochastic contribution, :math:`\text{prefactor}\sqrt{2T\boldsymbol{\mathcal{M}}}d\boldsymbol{W}`, where :math:`\boldsymbol{\mathcal{M}}` is the mobility matrix and :math:`d\boldsymbol{W}` is a Wiener process. +template auto call_mdot(Solver &myself, pyarray_c &forces) { + int N = myself.getNumberParticles(); + if (forces.size() < 3 * N and forces.size() > 0) { + throw std::runtime_error("The forces array must have size 3*N."); + } + auto f = forces.size() ? cast_to_const_real(forces) : nullptr; + auto result = + py::array_t(py::array::ShapeContainer({3 * N})); + result.attr("fill")(0); + myself.Mdot(f, cast_to_real(result)); + return result.reshape({N, 3}); +} + +const char *mdot_docstring = R"pbdoc( +Computes the product of the Mobility matrix with a group of forces, :math:`\boldsymbol{\mathcal{M}}\boldsymbol{F}`. It is required that :py:mod:`setPositions` has been called before calling this function. -Both inputs must have precision given by the precision attribute of the module. -Both inputs must have size 3*N, where N is the number of particles. -The arrays are ordered as :code:`[f0x, f0y, f0z, f1x, f1y, f1z, ...]`. Parameters ---------- -result : array_like, - Where the result will be stored. The result will have the same format as the forces array. -prefactor : float, optional - Prefactor to multiply the result by. Default is 1.0. +forces : array_like, + Forces acting on the particles. Must have shape (N, 3), where N is the number of particles. + +Returns +------- +array_like + The result of the product. The result will have the same format as the forces array. )pbdoc"; +template +void call_initialize(Solver &myself, libmobility::real T, libmobility::real eta, + libmobility::real a, int N, libmobility::real tol) { + libmobility::Parameters par; + par.temperature = T; + par.viscosity = eta; + par.hydrodynamicRadius = {a}; + par.tolerance = tol; + par.numberParticles = N; + myself.initialize(par); +} + +template void call_setPositions(Solver &myself, pyarray_c &pos) { + myself.setPositions(cast_to_const_real(pos)); +} +template +auto call_hydrodynamicVelocities(Solver &myself, pyarray_c &forces, + libmobility::real prefactor) { + int N = myself.getNumberParticles(); + if (forces.size() < 3 * N and forces.size() > 0) { + throw std::runtime_error("The forces array must have size 3*N."); + } + auto f = forces.size() ? cast_to_const_real(forces) : nullptr; + auto result = py::array_t({3 * N}); + result.attr("fill")(0); + myself.hydrodynamicVelocities(f, cast_to_real(result), prefactor); + return result.reshape({N, 3}); +} const char *hydrodynamicvelocities_docstring = R"pbdoc( Computes the hydrodynamic (deterministic and stochastic) velocities. @@ -136,58 +200,57 @@ Parameters ---------- forces : array_like, optional Forces acting on the particles. -result : array_like - Where the result will be stored. + prefactor : float, optional Prefactor to multiply the result by. Default is 1.0. -)pbdoc"; +Returns +------- +array_like + The resulting velocities. Shape is (N, 3), where N is the number of particles. +)pbdoc"; -#define xMOBILITY_PYTHONIFY(MODULENAME, EXTRACODE, documentation) \ - PYBIND11_MODULE(MODULENAME, m){ \ - using real = libmobility::real; \ - using Parameters = libmobility::Parameters; \ - using Configuration = libmobility::Configuration; \ - auto solver = py::class_(m, MOBILITYSTR(MODULENAME), documentation); \ - solver.def(py::init([](std::string perx, std::string pery, std::string perz){ \ - return std::unique_ptr(new MODULENAME(createConfiguration(perx, pery, perz))); }), \ - constructor_docstring, "periodicityX"_a, "periodicityY"_a, "periodicityZ"_a). \ - def("initialize", [](MODULENAME &myself, real T, real eta, real a, int N, real tol){ \ - Parameters par; \ - par.temperature = T; \ - par.viscosity = eta; \ - par.hydrodynamicRadius = {a}; \ - par.tolerance = tol; \ - par.numberParticles = N; \ - myself.initialize(par); \ - }, \ - initialize_docstring, \ - "temperature"_a, "viscosity"_a, \ - "hydrodynamicRadius"_a, \ - "numberParticles"_a, \ - "tolerance"_a = 1e-4). \ - def("setPositions", [](MODULENAME &myself, pyarray pos){myself.setPositions(cast_to_const_real(pos));}, \ - "The module will compute the mobility according to this set of positions.", \ - "positions"_a). \ - def("Mdot", [](MODULENAME &myself, pyarray forces, pyarray result){\ - auto f = forces.size()?cast_to_const_real(forces):nullptr; \ - myself.Mdot(f, cast_to_real(result));}, \ - mdot_docstring, \ - "forces"_a = pyarray(), "result"_a). \ - def("sqrtMdotW", [](MODULENAME &myself, pyarray result, libmobility::real prefactor){ \ - myself.sqrtMdotW(cast_to_real(result), prefactor);}, \ - sqrtMdotW_docstring, \ - "result"_a = pyarray(), "prefactor"_a = 1.0). \ - def("hydrodynamicVelocities", [](MODULENAME &myself, pyarray forces,\ - pyarray result, libmobility::real prefactor){ \ - auto f = forces.size()?cast_to_const_real(forces):nullptr; \ - myself.hydrodynamicVelocities(f, cast_to_real(result), prefactor);}, \ - hydrodynamicvelocities_docstring, \ - "forces"_a = pyarray(), "result"_a = pyarray(), "prefactor"_a = 1). \ - def("clean", &MODULENAME::clean, "Frees any memory allocated by the module."). \ - def_property_readonly_static("precision", [](py::object){return MODULENAME::precision;}, R"pbdoc(Compilation precision, a string holding either float or double.)pbdoc"); \ - EXTRACODE\ +template +auto call_construct(std::string perx, std::string pery, std::string perz) { + return std::unique_ptr( + new Solver(createConfiguration(perx, pery, perz))); } -#define MOBILITY_PYTHONIFY(MODULENAME, documentationPrelude) xMOBILITY_PYTHONIFY(MODULENAME,; ,documentationPrelude) -#define MOBILITY_PYTHONIFY_WITH_EXTRA_CODE(MODULENAME, EXTRA, documentationPrelude) xMOBILITY_PYTHONIFY(MODULENAME, EXTRA, documentationPrelude) + +#define xMOBILITY_PYTHONIFY(MODULENAME, EXTRACODE, documentation) \ + PYBIND11_MODULE(MODULENAME, m) { \ + using real = libmobility::real; \ + using Parameters = libmobility::Parameters; \ + using Configuration = libmobility::Configuration; \ + auto solver = \ + py::class_(m, MOBILITYSTR(MODULENAME), documentation); \ + solver \ + .def(py::init(&call_construct), constructor_docstring, \ + "periodicityX"_a, "periodicityY"_a, "periodicityZ"_a) \ + .def("initialize", call_initialize, initialize_docstring, \ + "temperature"_a, "viscosity"_a, "hydrodynamicRadius"_a, \ + "numberParticles"_a, "tolerance"_a = 1e-4) \ + .def("setPositions", call_setPositions, \ + "The module will compute the mobility according to this set of " \ + "positions.", \ + "positions"_a) \ + .def("Mdot", call_mdot, mdot_docstring, \ + "forces"_a = pyarray()) \ + .def("sqrtMdotW", call_sqrtMdotW, sqrtMdotW_docstring, \ + "prefactor"_a = 1.0) \ + .def("hydrodynamicVelocities", \ + call_hydrodynamicVelocities, \ + hydrodynamicvelocities_docstring, "forces"_a = pyarray_c(), \ + "prefactor"_a = 1) \ + .def("clean", &MODULENAME::clean, \ + "Frees any memory allocated by the module.") \ + .def_property_readonly_static( \ + "precision", [](py::object) { return MODULENAME::precision; }, \ + R"pbdoc(Compilation precision, a string holding either float or double.)pbdoc"); \ + EXTRACODE \ + } +#define MOBILITY_PYTHONIFY(MODULENAME, documentationPrelude) \ + xMOBILITY_PYTHONIFY(MODULENAME, ;, documentationPrelude) +#define MOBILITY_PYTHONIFY_WITH_EXTRA_CODE(MODULENAME, EXTRA, \ + documentationPrelude) \ + xMOBILITY_PYTHONIFY(MODULENAME, EXTRA, documentationPrelude) #endif diff --git a/solvers/SelfMobility/mobility.h b/solvers/SelfMobility/mobility.h index 96d2f74..d78bd0f 100644 --- a/solvers/SelfMobility/mobility.h +++ b/solvers/SelfMobility/mobility.h @@ -52,6 +52,9 @@ class SelfMobility: public libmobility::Mobility{ void setPositions(const real* ipositions) override{ } void Mdot(const real* forces, real* result) override{ + if(not forces){ + std::fill(result, result+3*numberParticles, 0); + } for(int i = 0; i<3*numberParticles; i++){ result[i] = forces[i]*selfMobility; } diff --git a/tests/test_fluctuation_dissipation.py b/tests/test_fluctuation_dissipation.py index e780e5d..82a9841 100644 --- a/tests/test_fluctuation_dissipation.py +++ b/tests/test_fluctuation_dissipation.py @@ -31,7 +31,7 @@ def fluctuation_dissipation_KS(M, fluctuation_method): # within +/- mu_a of the real mean # with probability mu_alpha for all N components N = Sigma.size - mu_alpha = 0.99 ** (1 / N) + mu_alpha = 0.999 ** (1 / N) mu_a = 0.05 Ns = int(round(2 * (norm.ppf(mu_alpha) / mu_a) ** 2)) ScaledNoise = np.full((N, Ns), np.nan) @@ -82,8 +82,6 @@ def test_fluctuation_dissipation( M = compute_M(solver, numberParticles) def fluctuation_method(): - sqrtmnoise = np.zeros(numberParticles * 3).astype(precision) - solver.sqrtMdotW(sqrtmnoise, prefactor=1.0) - return sqrtmnoise + return solver.sqrtMdotW(prefactor=1.0).flatten() fluctuation_dissipation_KS(M, fluctuation_method) diff --git a/tests/test_interface.py b/tests/test_interface.py index b50bbaf..2a3475e 100644 --- a/tests/test_interface.py +++ b/tests/test_interface.py @@ -1,16 +1,25 @@ import pytest from libMobility import * import numpy as np +from utils import sane_parameters @pytest.mark.parametrize( ("Solver", "periodicity"), - [(SelfMobility, ("open", "open", "open"))], + [ + (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")), + ], ) def test_contiguous(Solver, periodicity): hydrodynamicRadius = 1.0 solver = Solver(*periodicity) - solver.setParameters(parameter=1) + parameters = sane_parameters[Solver.__name__] + solver.setParameters(**parameters) numberParticles = 1 solver.initialize( temperature=1.0, @@ -24,9 +33,110 @@ def test_contiguous(Solver, periodicity): positions = np.random.rand(numberParticles, 3).astype(precision) solver.setPositions(positions) - # Set positions so that they are not contiguous matrix = np.random.rand(numberParticles * 3, numberParticles * 3).astype(precision) positions = matrix[:, 0] - with pytest.raises(RuntimeError): - solver.setPositions(positions) + solver.setPositions(positions) + + +@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")), + ], +) +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, + ) + + # Set precision to be the same as compiled precision + precision = np.float32 if Solver.precision == "float" else np.float64 + positions = np.random.rand(numberParticles, 3).astype(precision) + forces = np.random.rand(numberParticles, 3).astype(precision) + solver.setPositions(positions) + mf = solver.Mdot(forces) + assert mf.shape == (numberParticles, 3) + forces = forces.reshape(3 * numberParticles) + mf = solver.Mdot(forces) + assert mf.shape == (numberParticles, 3) + + +@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")), + ], +) +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, + ) + + # Set precision to be the same as compiled precision + precision = np.float32 if Solver.precision == "float" else np.float64 + positions = np.random.rand(numberParticles, 3).astype(precision) + solver.setPositions(positions) + sqrtmw = solver.sqrtMdotW() + assert sqrtmw.shape == (numberParticles, 3) + + +@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")), + ], +) +def test_returns_hydrodisp(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, + ) + + # Set precision to be the same as compiled precision + precision = np.float32 if Solver.precision == "float" else np.float64 + positions = np.random.rand(numberParticles, 3).astype(precision) + forces = np.random.rand(numberParticles, 3).astype(precision) + solver.setPositions(positions) + sqrtmw = solver.hydrodynamicVelocities() + assert sqrtmw.shape == (numberParticles, 3) + sqrtmw = solver.hydrodynamicVelocities(forces) + assert sqrtmw.shape == (numberParticles, 3) diff --git a/tests/test_mobility_matrix.py b/tests/test_mobility_matrix.py index 870f215..f8ddf2e 100644 --- a/tests/test_mobility_matrix.py +++ b/tests/test_mobility_matrix.py @@ -55,9 +55,8 @@ def test_self_mobility_selfmobility(): ) positions = np.zeros((1, 3), dtype=precision) solver.setPositions(positions) - result = np.zeros(3, dtype=precision) forces = np.ones(3, dtype=precision) - solver.Mdot(forces, result) + result = solver.Mdot(forces) m0 = 1.0 / (6 * np.pi * viscosity * hydrodynamicRadius) assert np.allclose(result, m0 * forces, rtol=0, atol=1e-7) @@ -79,9 +78,8 @@ def test_self_mobility_nbody(): ) positions = np.zeros((1, 3), dtype=precision) solver.setPositions(positions) - result = np.zeros(3, dtype=precision) forces = np.ones(3, dtype=precision) - solver.Mdot(forces, result) + result = solver.Mdot(forces) m0 = 1.0 / (6 * np.pi * viscosity * hydrodynamicRadius) assert np.allclose(result, m0 * forces, rtol=0, atol=1e-7) @@ -111,9 +109,8 @@ def test_self_mobility_pse_cubic_box(psi): ) positions = np.zeros((1, 3), dtype=precision) solver.setPositions(positions) - result = np.zeros(3, dtype=precision) forces = np.ones(3, dtype=precision) - solver.Mdot(forces, result) + result = solver.Mdot(forces) leff = hydrodynamicRadius / parameters["Lx"] m0 = ( 1.0 diff --git a/tests/test_precision.py b/tests/test_precision.py index 2256135..5670504 100644 --- a/tests/test_precision.py +++ b/tests/test_precision.py @@ -37,8 +37,7 @@ def test_precision(Solver, periodicity): solver.setPositions(positions) forces = np.ones(size, dtype=precision) - mf = np.zeros(size, dtype=precision) - solver.Mdot(forces, mf) + mf = solver.Mdot(forces) zeros = np.zeros(size, dtype=precision) assert ( @@ -51,7 +50,7 @@ def test_precision(Solver, periodicity): [(SelfMobility, ("open", "open", "open"))], ) def test_incorrect_precision(Solver, periodicity): - # tests that Mdot gives an error when inputs are in the wrong precision + # libMobility should work even if the inputs have an unexpected precision hydrodynamicRadius = 1.0 @@ -71,11 +70,8 @@ def test_incorrect_precision(Solver, periodicity): positions = np.random.rand(numberParticles, 3).astype(precision_bad) size = 3 * numberParticles - with pytest.raises(RuntimeError): - solver.setPositions(positions) + solver.setPositions(positions) positions = positions.astype(precision_good) solver.setPositions(positions) forces = np.ones(size, dtype=precision_bad) - mf = np.zeros(size, dtype=precision_bad) - with pytest.raises(RuntimeError): - solver.Mdot(forces, mf) + mf = solver.Mdot(forces) diff --git a/tests/utils.py b/tests/utils.py index 0f0ad3c..ac1da48 100644 --- a/tests/utils.py +++ b/tests/utils.py @@ -16,9 +16,7 @@ def compute_M(solver, numberParticles): I = np.identity(size, dtype=precision) for i in range(0, size): forces = I[:, i].copy() - mf = np.zeros(size, dtype=precision) - solver.Mdot(forces, mf) - M[:, i] = mf + M[:, i] = solver.Mdot(forces).reshape(3 * numberParticles) return M