Adding tests

.github/workflows/ci.yml

@@ -57,6 +57,13 @@ jobs:
         python -c "import libMobility"
         python -c "from libMobility import DPStokes, PSE, NBody, SelfMobility"
 
+    - name: Run pytest
+      run: |
+        micromamba activate libmobility
+        pytest -sv tests/test_initialization.py tests/test_precision.py
+        #Only SelfMobility tests can run in the CPU
+        pytest -sv -k SelfMobility tests/test*py
+
     - name: Install docs dependencies
       run: |
         pip install -r docs/requirements.txt

environment.yml

@@ -12,3 +12,5 @@ dependencies:
   - pybind11
   - numpy
   - python 3.11.*
+  - pytest
+  - scipy

include/MobilityInterface/MobilityInterface.h

@@ -106,7 +106,6 @@ namespace libmobility{
     }
 
     //Equivalent to calling Mdot and then stochasticDisplacements, can be faster in some solvers
-    // Donev: The README says there is a torque parameter here but I do not see it (I asked that it be removed but maybe this README never got updated)?
     virtual void hydrodynamicVelocities(const real* forces, real* result, real prefactor = 1){
       Mdot(forces, result);
       sqrtMdotW(result, prefactor);

include/MobilityInterface/pythonify.h

@@ -8,7 +8,7 @@ The MOBILITY_PYTHONIFY(className, description) macro creates a pybind11 module f
 #include<pybind11/numpy.h>
 namespace py = pybind11;
 using namespace pybind11::literals;
-using pyarray = py::array_t<libmobility::real>;
+using pyarray = py::array;
 
 #define MOBILITYSTR(s) xMOBILITYSTR(s)
 #define xMOBILITYSTR(s) #s
@@ -31,6 +31,26 @@ inline auto createConfiguration(std::string perx, std::string pery, std::string
   return conf;
 }
 
+template<typename T>
+void check_dtype(pyarray &arr){
+  if(not py::isinstance<py::array_t<T>>(arr)){
+    throw std::runtime_error("Input array must have the correct data type.");
+  }
+  if(not py::isinstance<py::array_t<T, py::array::c_style | py::array::forcecast>>(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){
+  check_dtype<libmobility::real>(arr);
+  return static_cast<libmobility::real*>(arr.mutable_data());
+}
+
+const libmobility::real* cast_to_const_real(pyarray &arr){
+  check_dtype<const libmobility::real>(arr);
+  return static_cast<const libmobility::real*>(arr.data());
+}
 const char *constructor_docstring = R"pbdoc(
 Initialize the module with a given set of periodicity conditions.
 
@@ -64,6 +84,8 @@ hydrodynamicRadius : float
 		Hydrodynamic radius of the particles.
 numberParticles : int
 		Number of particles in the system.
+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(
@@ -130,35 +152,37 @@ prefactor : float, optional
   solver.def(py::init([](std::string perx, std::string pery, std::string perz){	\
     return std::unique_ptr<MODULENAME>(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){ \
+    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).						\
-    def("setPositions", [](MODULENAME &myself, pyarray pos){myself.setPositions(pos.data());}, \
+	  "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()?forces.data():nullptr;			\
-      myself.Mdot(f, result.mutable_data());},			\
-      mdot_docstring,	\
+      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(result.mutable_data(), prefactor);}, \
-      sqrtMdotW_docstring, \
+      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()?forces.data():nullptr;			\
-      myself.hydrodynamicVelocities(f, result.mutable_data(), prefactor);}, \
-	hydrodynamicvelocities_docstring, \
+      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"); \

solvers/DPStokes/mobility.h

@@ -23,6 +23,7 @@ class DPStokes: public libmobility::Mobility{
   real lanczosTolerance;
   std::uint64_t lanczosSeed;
   std::shared_ptr<LanczosStochasticVelocities> lanczos;
+  std::string wallmode;
 public:
 
   DPStokes(Configuration conf){
@@ -33,31 +34,13 @@ class DPStokes: public libmobility::Mobility{
 	    conf.periodicityZ == libmobility::periodicity_mode::two_walls)
 	    )
       throw std::runtime_error("[DPStokes] This is a doubly periodic solver");
+    if(conf.periodicityZ == periodicity_mode::open) wallmode = "nowall";
+    else if(conf.periodicityZ == periodicity_mode::single_wall) wallmode = "bottom";
+    else if(conf.periodicityZ == periodicity_mode::two_walls) wallmode = "slit";
   }
 
   void setParametersDPStokes(DPStokesParameters i_dppar){
     this->dppar = i_dppar;
-    // DPStokesParameters dppar;
-    // int nx = -1;
-    // int ny = -1;
-    // int nz = -1;
-    // dppar.dt
-    // dppar.viscosity
-    // dppar.Lx
-    // dppar.Ly
-    // dppar.zmin
-    //   dppar. zmax
-    // //Tolerance will be ignored in DP mode, TP will use only tolerance and nxy/nz
-    // dppar.tolerance = 1e-5;
-    // dppar.w
-    // dppar.w_d;
-    // dppar.hydrodynamicRadius = -1;
-    // dppar.beta = -1;
-    // dppar.beta_d = -1;
-    // dppar.alpha = -1;
-    // dppar.alpha_d = -1;
-    // //Can be either none, bottom, slit or periodic
-    // dppar.mode;
     dpstokes = std::make_shared<uammd_dpstokes::DPStokesGlue>();
   }
 
@@ -66,6 +49,25 @@ class DPStokes: public libmobility::Mobility{
     this->dppar.viscosity = ipar.viscosity;
     this->temperature = ipar.temperature;
     this->lanczosTolerance = ipar.tolerance;
+    this->dppar.mode = this->wallmode;
+    this->dppar.hydrodynamicRadius = ipar.hydrodynamicRadius[0];
+    this->dppar.w = 6;
+    this->dppar.beta = 1.714*this->dppar.w;
+    real h = this->dppar.hydrodynamicRadius/1.554;
+    this->dppar.alpha = this->dppar.w/2.0;
+    this->dppar.tolerance = 1e-6;
+    this->dppar.nx = int(this->dppar.Lx/h + 0.5);
+    this->dppar.ny = int(this->dppar.Ly/h + 0.5);
+    // Add a buffer of w*h/2 when there is an open boundary
+    if(this->wallmode == "nowall"){
+      this->dppar.zmax += this->dppar.w*h/2;
+      this->dppar.zmin -= this->dppar.w*h/2;
+    }
+    if(this->wallmode == "bottom"){
+      this->dppar.zmin -= this->dppar.w*h/2;
+    }
+    real Lz = this->dppar.zmax - this->dppar.zmin;
+    this->dppar.nz = M_PI*Lz/(h);
     dpstokes->initialize(dppar, this->numberParticles);
     Mobility::initialize(ipar);
   }
@@ -78,22 +80,6 @@ class DPStokes: public libmobility::Mobility{
     dpstokes->Mdot(forces, nullptr, result, nullptr);
   }
 
-  void sqrtMdotW(real* result, real prefactor = 1) override{
-    if(this->temperature == 0) return;
-    if(not dpstokes)
-      throw std::runtime_error("[libMobility] You must initialize the base class in order to use the default stochastic displacement computation");
-    if(not lanczos){
-      if(this->lanczosSeed==0){//If a seed is not provided, get one from random device
-	this->lanczosSeed = std::random_device()();
-      }
-      lanczos = std::make_shared<LanczosStochasticVelocities>(this->numberParticles,
-							      this->lanczosTolerance, this->lanczosSeed);
-    }
-    lanczos->sqrtMdotW([this](const real*f, real* mv){Mdot(f, mv);}, result, prefactor);
-  }
-
-  // Donev: Is there a Lanczos implementation on the GPU or that is always done on CPU? If there is a GPU version, is it possible with the current interface to make a more efficient hydrodynamicDisplacements routine that does not do the CPU<->GPU twice of say the positions. That is, there is no return back to the GPU until both Mdot and sqrtMdotW are finished? I cannot follow the UAMMD stuff so just tell me what the code does now and if it is possible to optimize further.
-
   void clean() override{
     Mobility::clean();
     dpstokes->clear();

solvers/DPStokes/python_wrapper.cpp

@@ -6,12 +6,12 @@ using DPStokesParameters = uammd_dpstokes::PyParameters;
 static const char *docstring = R"pbdoc(
 In the Doubly periodic Stokes geometry (DPStokes), an incompressible fluid exists in a domain which is periodic in the plane and open (or walled) in the third direction.
 
+When the periodicity is set to :code:`single_wall` a wall in the bottom of the domain is added.
+
 Parameters
 ----------
 dt : float
 		The time step.
-viscosity : float
-		The viscosity of the fluid.
 Lx : float
 		The box size in the x direction.
 Ly : float
@@ -20,49 +20,21 @@ zmin : float
 		The minimum value of the z coordinate.
 zmax : float
 		The maximum value of the z coordinate.
-w : float
-		The support of the particle force spreading kernel.
-w_d : float
-		The support of the dipole force spreading kernel.
-hydrodynamicRadius : float
-		The hydrodynamic radius of the particle.
-beta : float
-		The width of the Exponential of the Semicircle (ES) particle force spreading kernel.
-beta_d : float
-		The width of the Exponential of the Semicircle (ES) particle dipole spreading kernel.
-alpha : float
-		Distance normalization factor for the particle force spreading kernel.
-alpha_d : float
-		Distance normalization factor for the dipole force spreading kernel.
-mode : str
-		The wall mode. Options are 'slit', 'bottom' and 'nowall'.
 )pbdoc";
 
-MOBILITY_PYTHONIFY_WITH_EXTRA_CODE(
-    DPStokes,
-    solver.def(
-        "setParameters",
-        [](DPStokes &self, real dt, real viscosity, real Lx, real Ly, real zmin,
-           real zmax, real w, real w_d, real hydrodynamicRadius, real beta,
-           real beta_d, real alpha, real alpha_d, std::string mode) {
-	  DPStokesParameters params;
-	  params.dt = dt;
-	  params.viscosity = viscosity;
-	  params.Lx = Lx;
-	  params.Ly = Ly;
-	  params.zmin = zmin;
-	  params.zmax = zmax;
-	  params.w = w;
-	  params.w_d = w_d;
-	  params.hydrodynamicRadius = hydrodynamicRadius;
-	  params.beta = beta;
-	  params.beta_d = beta_d;
-	  params.alpha = alpha;
-	  params.alpha_d = alpha_d;
-	  params.mode = mode;
-          self.setParametersDPStokes(params);
-        },
-	"dt"_a, "viscosity"_a, "Lx"_a, "Ly"_a, "zmin"_a, "zmax"_a, "w"_a,
-	"w_d"_a, "hydrodynamicRadius"_a, "beta"_a, "beta_d"_a, "alpha"_a,
-	"alpha_d"_a, "mode"_a);,
-    docstring);
+MOBILITY_PYTHONIFY_WITH_EXTRA_CODE(DPStokes,
+                                   solver.def(
+                                       "setParameters",
+                                       [](DPStokes &self, real dt, real Lx,
+                                          real Ly, real zmin, real zmax) {
+                                         DPStokesParameters params;
+                                         params.dt = dt;
+                                         params.Lx = Lx;
+                                         params.Ly = Ly;
+                                         params.zmin = zmin;
+                                         params.zmax = zmax;
+                                         self.setParametersDPStokes(params);
+                                       },
+                                       "dt"_a, "Lx"_a, "Ly"_a, "zmin"_a,
+                                       "zmax"_a);
+                                   , docstring);