Add torques (#21)

* Add torques to pythonify * Add torque to mobility interface * Add torques to DPStokes * Add torques to NBody * Add torques to PSE * Add torques to selfMobility * Update pythonify * Update python example * Update cpp example * Update test * Update test * Update test * Update test * Update test * Typo * Fix initialize * Typo * Update pythonify * Update pythonify * Typo * Update sqrtM * Update * Update * Typo * Change name * Update * Always send valid pointers to hydrodynamicVelocities * dont complain if receiving angular * Typo in tests * added torques to selfmobility test * angular mobility matrix test setup * added torques to selfmobility * updated interface so that the linear fluctuation dissipation test is passing. the problem was that a zero vector of angular velocities was always being passed into sqrtMdotW inside pythonify.h so solvers where torques weren't implemented would always throw the error, even though we were only trying to compute with forces. * removed unnecessary allocation for when torques aren't being used in sqrtMdotW * Squashed commit of the following: commit ac2b405 Author: Ryker Fish <rfish@mines.edu> Date: Tue Jul 16 19:41:16 2024 -0600 regenerated data w new defaults after verifying against old repository commit 07b4794 Author: Ryker Fish <rfish@mines.edu> Date: Tue Jul 16 19:05:08 2024 -0600 added parameter to initialize that sets whether the box size can get altered or not commit 447e543 Author: Ryker Fish <rfish@mines.edu> Date: Tue Jul 16 18:22:27 2024 -0600 added mode that preserves L commit 74ec436 Author: Ryker Fish <rfish@mines.edu> Date: Tue Jul 16 17:23:24 2024 -0600 checkpoint commit- adding functionality to have parameter select that modifies the box size and another version that preserves it commit 1113390 Author: Ryker Fish <rfish@mines.edu> Date: Wed Jul 10 14:59:33 2024 -0600 the changes to DPStokes parameters seem to made the mobility matrix less symmetric commit 800f634 Author: Ryker Fish <rfish@mines.edu> Date: Wed Jul 10 14:56:08 2024 -0600 had to change the ref file to use the right positions. nbody test working. commit 830866e Author: Ryker Fish <rfish@mines.edu> Date: Mon Jul 8 18:46:51 2024 -0600 the two pair DPStokes tests are passing. currently having issues with the nbody test commit 9b3a4d3 Author: Ryker Fish <rfish@mines.edu> Date: Fri Jun 28 14:19:28 2024 -0600 self mobility tests done commit c4f2664 Author: Ryker Fish <rfish@mines.edu> Date: Wed Jun 12 15:31:55 2024 -0600 larger safety factor near open boundary for dpstokes commit 5da7450 Author: Ryker Fish <rfish@mines.edu> Date: Wed Jun 12 12:15:33 2024 -0600 modify how N gets rounded in all dimensions to be fft friendly commit 51db36a Author: Ryker Fish <rfish@mines.edu> Date: Wed Jun 12 11:58:39 2024 -0600 changed selection of grid parameters to not modify h. instead, modify the box size to make L/h an integer commit 5127394 Author: Ryker Fish <rfish@mines.edu> Date: Thu Jun 6 14:07:01 2024 -0600 changed test pass threshold commit 4fb82bf Author: Ryker Fish <rfish@mines.edu> Date: Thu Jun 6 13:42:52 2024 -0600 regenerated all ref data using libmobility implementations in single precision commit 4392903 Author: Ryker Fish <rfish@mines.edu> Date: Thu Jun 6 13:25:04 2024 -0600 regenerated refs using gpu implementations from dpstokes library. tests are now passing if the same kernel parameters are used as the ref implementation commit 1ce606a Author: Ryker Fish <rfish@mines.edu> Date: Thu May 30 16:06:32 2024 -0600 fixed bug in nbody wall correction that came from normalizing by a factor of 1/8pi in the original implementation instead of 1/6pi in the Swan paper. commit edf78f5 Author: Ryker Fish <rfish@mines.edu> Date: Thu May 30 16:05:21 2024 -0600 switched all tests to w=4 kernels, added nbody to pair mobility tests commit c575ef2 Author: Ryker Fish <rfish@mines.edu> Date: Wed May 29 15:11:48 2024 -0600 change dpstokes to use w=4 kernel by default for forces commit 070728b Author: Ryker Fish <rfish@mines.edu> Date: Tue May 28 16:57:09 2024 -0600 started pair mobility test commit ca4e1c7 Author: Ryker Fish <rfish@mines.edu> Date: Sat May 18 12:19:26 2024 -0600 combined self mobility tests into one function commit e90e9e2 Author: Ryker Fish <rfish@mines.edu> Date: Fri May 17 20:45:24 2024 -0600 bug fix in NBody with wall correction- test now passing. the fix uses the correct constants from the Swan correction paper- the previous implementation was using constants based on a self mobility constant of 1/(8*pi*eta*a) instead of 1/(6*pi*eta*a) commit 98e2ae0 Author: Ryker Fish <rfish@mines.edu> Date: Fri May 17 15:37:25 2024 -0600 changed dpstokes test to compare to dpstokes reference and added a not-working nbody test commit 7aaa0b1 Merge: 79fe6a7 d6144a5 Author: Ryker <rfish@mines.edu> Date: Mon May 13 14:25:35 2024 -0600 Merge pull request #18 from stochasticHydroTools/minor_test_edits Minor test edits commit d6144a5 Author: Ryker Fish <rfish@mines.edu> Date: Thu May 9 15:40:00 2024 -0600 removed test cases for small numbers of particles for the fluctuation dissipation test. commit 0e368d4 Author: Ryker <rfish@mines.edu> Date: Thu May 9 14:37:44 2024 -0600 Edit comment on nbody test Comment had the answer for the PSE test. Removed the periodic correction part commit 79fe6a7 Author: Ryker Fish <rfish@mines.edu> Date: Thu May 9 12:47:42 2024 -0600 now testing if mobility is 0 on the bottom wall commit 3c5002e Author: Ryker Fish <rfish@mines.edu> Date: Wed May 8 19:02:35 2024 -0600 possible correction to DPStokes code to add a buffer the top instead of the bottom when there is a wall commit c66303c Author: Ryker Fish <rfish@mines.edu> Date: Wed May 8 18:27:40 2024 -0600 selfmobility test working in comparison to refernce results for one and two wall cases commit 70aa329 Author: Ryker Fish <rfish@mines.edu> Date: Wed May 8 17:12:19 2024 -0600 mobility doesn't go to zero at the wall unless this section is commented out. commit fe2c964 Author: Ryker Fish <rfish@mines.edu> Date: Wed May 8 16:37:04 2024 -0600 first pass at computing translational self mobility components for bottom wall scenario * modified wall tests to still pass in single precision by using variable tolerance based on machine epsilon * modified how symmetric matrix has to be for fluctuation dissipation tests to accomodate dpstokes parameters * changed to only allocate angular vec if torques are provided so that PSE and NBody work when torques aren't provided * interface test for mdot with torques * updated all interfaces to accept torques * started implementing NBody torques. naive algorithm works for main diagonal blocks- still need to do off diagonal blocks * fast algorithm now working for main diagonal * updated angular test to use different torques and forces * added off-diagonal blocks for torques and fixed bug in memory allocation for fast algorithm when using torques * added torque test for when forces are zero with nonzero torques * added test for freespace nbody with torques * bug fix in nbody UT component * add test for PSE not supporting torques * added torque parameters * add DPStokes to mob. matrix tests and remove PSE from tests involving torques * self mobility tests for DPstokes with torques * dpstokes pair wall tests complete * rewrite linear wall test to be more like angular test * fix bug in NBody bloc kalgorithm * added parameterization to test different nbody algorithms to show bug * bug fix- thanks Raul! * add block algorithm back to nbody tests * added torques to nbody block * added nbody wall kernels for self mobilities. still need to do pair mobilities * checkpoint commit for today, stuff still kinda broken but I think I know why * simplified tests since ref files have changed * added wall test for particle offset in y direction * wall kernels for torques in Nbody now working after bug fixes * added some comments on the subtle parts of the NBody wall kernels * minor optimization changes (went division hunting) and standardized notation between functions * add line to example to indicate to compile in double if necessary * change python example function calls to reflect interface * Set pybind11 as required * Take some code out of the PYTHONIFY macro * Format * Fix module name * Fix typo in docstring * Set C++ std to 17 * Update docstrings * Check if torques are being passed when needsTorques=False * Fix test not setting needsTorque=True * Install test scripts to /share/ * Add a different test per file * Improve mkl detection --------- Co-authored-by: Ryker Fish <rfish@mines.edu>
stochasticHydroTools · Sep 19, 2024 · f07856c · f07856c
1 parent d013acf
commit f07856c
Show file tree

Hide file tree

Showing 32 changed files with 1,044 additions and 246 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -5,7 +5,7 @@ set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/" "$ENV{CO
 set(CMAKE_BUILD_TYPE Release)
 add_compile_definitions(PUBLIC MAXLOGLEVEL=3)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
-set(CMAKE_CXX_STANDARD 14)
+set(CMAKE_CXX_STANDARD 17)
 set(CMAKE_CUDA_STANDARD 14)
 set(CMAKE_CUDA_STANDARD_REQUIRED ON)
 set(CMAKE_CUDA_SEPARABLE_COMPILATION OFF)
@@ -25,8 +25,7 @@ find_package(LAPACK)
 if(BLAS_FOUND AND LAPACK_FOUND)
   message(STATUS "BLAS and LAPACK found")
   # Additional check for MKL
-  find_path(MKL_INCLUDE_DIR NAMES mkl.h PATHS $ENV{MKLROOT}/include)
-  if(MKL_INCLUDE_DIR)
+  if(LAPACK_LIBRARIES MATCHES "mkl")
     message(STATUS "Intel MKL environment detected")
     add_compile_definitions(PUBLIC USE_MKL)
     include_directories(${MKL_INCLUDE_DIR})
@@ -58,7 +57,7 @@ else()
 endif()
 
 find_package(Python3 REQUIRED COMPONENTS Interpreter Development)
-find_package(pybind11 CONFIG)
+find_package(pybind11 REQUIRED CONFIG)
 include_directories(${pybind11_INCLUDE_DIRS})
 include_directories(${Python3_INCLUDE_DIRS})
 include_directories(${UAMMD_INCLUDE} ${BLAS_INCLUDE_DIRS} ${LAPACKE_INCLUDE_DIRS})
@@ -89,3 +88,13 @@ ENDIF()
 install(DIRECTORY include/MobilityInterface DESTINATION ${CMAKE_INSTALL_PREFIX}/include/)
 install(DIRECTORY include/third_party/LanczosAlgorithm/include/ DESTINATION ${CMAKE_INSTALL_PREFIX}/include/MobilityInterface/)
 install(DIRECTORY include/third_party/uammd/src/ DESTINATION ${CMAKE_INSTALL_PREFIX}/include/uammd/src)
+
+
+enable_testing()
+install(DIRECTORY tests/ DESTINATION ${CMAKE_INSTALL_PREFIX}/share/libMobility/tests/ FILES_MATCHING PATTERN "*.py" PATTERN "*.npz")
+# Add a test for each Python file in the tests directory
+file(GLOB TESTS tests/*.py)
+foreach(TEST ${TESTS})
+	get_filename_component(TEST_NAME ${TEST} NAME_WE)
+	add_test(NAME ${TEST_NAME} COMMAND pytest -vs ${TEST})
+endforeach()
diff --git a/examples/cpp/CMakeLists.txt b/examples/cpp/CMakeLists.txt
@@ -9,6 +9,8 @@ set(CMAKE_CUDA_STANDARD 14)
 set(CMAKE_CUDA_STANDARD_REQUIRED ON)
 set(CMAKE_CUDA_SEPARABLE_COMPILATION OFF)
 set(UAMMD_INCLUDE include/third_party/uammd/src include/third_party/uammd/src/third_party)
+# enable the line below if libmobility is compiled in double
+# add_compile_definitions(PUBLIC DOUBLE_PRECISION)
 set(BLA_STATIC OFF)
 set(BLA_VENDOR Intel10_64lp)
 find_package(BLAS)

diff --git a/examples/cpp/example.cpp b/examples/cpp/example.cpp
@@ -51,7 +51,7 @@ auto computeMFWithSolver(std::shared_ptr<MobilityBase> solver,
 			 std::vector<scalar> &forces){
   std::vector<scalar> result(pos.size(), 0);
   solver->setPositions(pos.data());
-  solver->Mdot(forces.data(), result.data());
+  solver->Mdot(forces.data(), nullptr, result.data(), nullptr);
   return result;
 }
 
@@ -74,6 +74,7 @@ int main(){
   par.numberParticles = numberParticles;
   par.tolerance = 1e-4;
   par.temperature = 1.0;
+  par.needsTorque = false;
 
   //Create two different solvers
   auto solver_pse = initializeSolver<PSE>(par);
@@ -86,7 +87,7 @@ int main(){
   //The solvers can be used to compute stochastic displacements, even if they do not provide a specific way to compute them (defaults to using the lanczos algorithm
   std::vector<scalar> noiseNBody(pos.size(), 0);
   scalar prefactor = 1.0;
-  solver_nbody->sqrtMdotW(noiseNBody.data(), prefactor);
+  solver_nbody->sqrtMdotW(noiseNBody.data(), nullptr, prefactor);
 
   //Remember to clean up when done
   solver_nbody->clean();

diff --git a/examples/python/example.py b/examples/python/example.py
@@ -1,13 +1,15 @@
-'''
+"""
 Raul P. Pelaez 2021. libMobility python interface usage example.
 One of the available solvers is chosen and then, using the same interface the solver is initialized and the mobility is applied.
 
-'''
+"""
+
 import numpy as np
 import libMobility
+
 # can also import modules individually by name, e.g.
 # from libMobility import SelfMobility
-# from libMobility import PSE 
+# from libMobility import PSE
 # ... etc
 
 # Each solver has its own help available, and libMobility itself offers information about the common interface.
@@ -18,35 +20,40 @@
 
 
 numberParticles = 3
-#Each module can be compiled using a different precision, you can know which one with this
-precision = np.float32 if libMobility.SelfMobility.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)
+# Each module can be compiled using a different precision, you can know which one with this
+precision = np.float32 if libMobility.SelfMobility.precision == "float" else np.float64
+pos = np.random.rand(3 * numberParticles).astype(precision)
+force = np.ones(3 * numberParticles).astype(precision)
 
-#Any solver can be used interchangeably, some of them need additional initialization via a "setParameters" function
+# Any solver can be used interchangeably, some of them need additional initialization via a "setParameters" function
 
-nb = libMobility.SelfMobility(periodicityX='open',periodicityY='open',periodicityZ='open')
+nb = libMobility.SelfMobility(
+    periodicityX="open", periodicityY="open", periodicityZ="open"
+)
 # to call with the alternate import, use the below
 # nb = SelfMobility(periodicityX='open',periodicityY='open',periodicityZ='open')
 
 # For NBody periodicityZ can also be single_wall
 # nb = libMobility.NBody(periodicityX='open',periodicityY='open',periodicityZ='open')
-# nb.setParametersNBody(algorithm="advise", Nbatch=1, NperBatch=numberParticles)
+# nb.setParameters(algorithm="advise", Nbatch=1, NperBatch=numberParticles)
 
 # nb = libMobility.PSE(periodicityX='periodic',periodicityY='periodic',periodicityZ='periodic')
-# nb.setParametersPSE(psi=1,   Lx=128, Ly=128, Lz=128,shearStrain=1)
-
-nb.initialize(temperature=1.0, viscosity = 1/(6*np.pi),
-              hydrodynamicRadius = 1.0,
-              numberParticles = numberParticles)
+# nb.setParameters(psi=1,   Lx=128, Ly=128, Lz=128,shearStrain=1)
+
+nb.initialize(
+    temperature=1.0,
+    viscosity=1 / (6 * np.pi),
+    hydrodynamicRadius=1.0,
+    numberParticles=numberParticles,
+    needsTorque=False,
+)
 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 = prefactor*sqrt(2*temperature*M)*dW
+result = nb.sqrtMdotW(prefactor=1.0)
 print("sqrt(2*T*M)*N(0,1):", result)
 nb.clean()
diff --git a/include/MobilityInterface/MobilityInterface.h b/include/MobilityInterface/MobilityInterface.h
@@ -22,6 +22,7 @@ namespace libmobility{
     real tolerance = 1e-4; //Tolerance for Lanczos fluctuations
     int numberParticles = -1;
     std::uint64_t seed = 0;
+    bool needsTorque = false;
   };
 
   //A list of parameters that cannot be changed by reinitializing a solver and/or are properties of the solver.
@@ -41,7 +42,9 @@ namespace libmobility{
     std::uint64_t lanczosSeed;
     real lanczosTolerance;
     std::shared_ptr<LanczosStochasticVelocities> lanczos;
+    std::vector<real> lanczosOutput;
     real temperature;
+    bool needsTorque = false;
   protected:
     Mobility(){};
   public:
@@ -81,45 +84,64 @@ namespace libmobility{
       this->lanczosSeed = par.seed;
       this->lanczosTolerance = par.tolerance;
       this->temperature = par.temperature;
+      this->needsTorque = par.needsTorque;
     }
 
     //Set the positions to construct the mobility operator from
     virtual void setPositions(const real* positions) = 0;
 
-    //Apply the mobility operator (M) to a series of forces (F) returns M*F
-    virtual void Mdot(const real* forces, real* result) = 0;
+    //Apply the grand mobility operator (\Omega) to a series of forces (F) and torques (T) to get the resulting linear (V) and angular (W) velocities
+    virtual void Mdot(const real* forces, const real* torques, real* linear, real* angular) = 0;
 
-    //Compute the stochastic displacements as result=prefactor*sqrt(M)*dW. Where dW is a vector of Gaussian random numbers
+    //Compute the stochastic displacements as result=prefactor*sqrt(\Omega)*dW. Where dW is a vector of Gaussian random numbers
     //If the solver does not provide a stochastic displacement implementation, the Lanczos algorithm will be used automatically
-    virtual void sqrtMdotW(real* result, real prefactor = 1){
+    virtual void sqrtMdotW(real* linear, real* angular, real prefactor = 1){
       if(this->temperature == 0) return;
+      if(prefactor == 0) return;
       if(not this->initialized)
 	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,
+	int numberElements = this->needsTorque ? (2*this->numberParticles) : this->numberParticles;
+	lanczos = std::make_shared<LanczosStochasticVelocities>(numberElements,
 								this->lanczosTolerance, this->lanczosSeed);
+	lanczosOutput.resize(3*numberElements);
       }
-      lanczos->sqrtMdotW([this](const real*f, real* mv){Mdot(f, mv);}, result, prefactor);
+      if(this->needsTorque and not angular)
+	throw std::runtime_error("[libMobility] This solver requires angular velocities when configured with torques");
+      std::fill(lanczosOutput.begin(), lanczosOutput.end(), 0);
+      lanczos->sqrtMdotW([this](const real*f, real* mv){
+	const real* t = this->needsTorque ? (f+3*this->numberParticles) : nullptr;
+        real* mt = this->needsTorque ? (mv+3*this->numberParticles) : nullptr;
+	Mdot(f, t, mv, mt);
+      }, lanczosOutput.data(), prefactor);
+      std::copy(lanczosOutput.begin(), lanczosOutput.begin()+3*this->numberParticles, linear);
+      if(this->needsTorque)
+	std::copy(lanczosOutput.begin()+3*this->numberParticles, lanczosOutput.end(), angular);
     }
 
     //Equivalent to calling Mdot and then stochasticDisplacements, can be faster in some solvers
-    virtual void hydrodynamicVelocities(const real* forces, real* result, real prefactor = 1){
-      if(forces){
-	Mdot(forces, result);
+    virtual void hydrodynamicVelocities(const real* forces, const real* torques, real* linear, real* angular, real prefactor = 1){
+      if(forces or torques){
+	Mdot(forces, torques, linear, angular);
       }
-      sqrtMdotW(result, prefactor);
+      sqrtMdotW(linear, angular, prefactor);
     }
 
     int getNumberParticles() const{
       return this->numberParticles;
     }
 
+    bool getNeedsTorque() const{
+      return this->needsTorque;
+    }
+
     //Clean any memory allocated by the solver
     virtual void clean(){
       lanczos.reset();
+      lanczosOutput.clear();
     }
   };
 }