espressomd · biermanncarl · Dec 16, 2018 · Dec 22, 2018 · Jan 1, 2019 · Jan 12, 2019
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -253,6 +253,23 @@ if(WITH_GSL)
   endif(GSL_FOUND)
 endif(WITH_GSL)
 
+find_package(BLAS)
+if(BLAS_FOUND)
+  set(BLAS 1)
+endif()
+find_package(LAPACK)
+if(LAPACK_FOUND)
+  set(LAPACK 1)
+endif()
+
+if(BLAS AND LAPACK)
+  set(STOKESIAN_DYNAMICS 1)
+endif()
+
+if(WITH_CUDA)
+  set(STOKESIAN_DYNAMICS_GPU 1)
+endif()
+
 if(WITH_VALGRIND_INSTRUMENTATION)
   find_package(PkgConfig)
   pkg_check_modules(VALGRIND valgrind)

diff --git a/cmake/FindCUDACompilerClang.cmake b/cmake/FindCUDACompilerClang.cmake
@@ -73,6 +73,14 @@ find_library(
   CUFFT_LIBRARY NAMES cufft PATHS ${CUDA_DIR}/lib64 ${CUDA_DIR}/lib
                                   /usr/local/nvidia/lib
                                   /usr/lib/x86_64-linux-gnu NO_DEFAULT_PATH)
+find_library(
+  CUDA_CUBLAS_LIBRARIES NAMES cublas PATHS ${CUDA_DIR}/lib64 ${CUDA_DIR}/lib
+                                    /usr/local/nvidia/lib
+                                    /usr/lib/x86_64-linux-gnu NO_DEFAULT_PATH)
+find_library(
+  CUDA_cusolver_LIBRARY NAMES cusolver PATHS ${CUDA_DIR}/lib64 ${CUDA_DIR}/lib
+                                  /usr/local/nvidia/lib
+                                  /usr/lib/x86_64-linux-gnu NO_DEFAULT_PATH)
 
 function(add_gpu_library)
   set(options STATIC SHARED MODULE EXCLUDE_FROM_ALL)

diff --git a/cmake/FindCUDACompilerHIP.cmake b/cmake/FindCUDACompilerHIP.cmake
@@ -53,6 +53,12 @@ list(APPEND HIP_HIPCC_FLAGS_RELWITHASSERT -O3 -g)
 
 find_library(ROCFFT_LIB name "rocfft" PATHS ${ROCM_HOME}/lib)
 
+find_library(ROCBLAS_LIB name "rocblas" PATHS ${ROCM_HOME}/lib)
+set(CUDA_CUBLAS_LIBRARIES ${ROCBLAS_LIB})
+
+find_library(ROCSOLVER_LIB name "rocsolver" PATHS ${ROCM_HOME}/lib)
+set(CUDA_cusolver_LIBRARY ${ROCSOLVER_LIB})
+
 function(add_gpu_library)
   hip_add_library(${ARGV})
   set(GPU_TARGET_NAME ${ARGV0})

diff --git a/cmake/cmake_config.cmakein b/cmake/cmake_config.cmakein
@@ -13,6 +13,10 @@
 
 #cmakedefine GSL
 
+#cmakedefine BLAS
+
+#cmakedefine LAPACK
+
 #cmakedefine VALGRIND_INSTRUMENTATION
 
 #define PACKAGE_NAME "${PROJECT_NAME}"

diff --git a/doc/doxygen/Doxyfile.in b/doc/doxygen/Doxyfile.in
@@ -790,6 +790,7 @@ WARN_LOGFILE           =
 # Note: If this tag is empty the current directory is searched.
 
 INPUT                  = "@CMAKE_SOURCE_DIR@/src" \
+                         "@CMAKE_SOURCE_DIR@/libs/stokesian_dynamics" \
                          "@CMAKE_SOURCE_DIR@/doc/doxygen/main.dox"
 
 # This tag can be used to specify the character encoding of the source files

diff --git a/doc/doxygen/bibliography.bib b/doc/doxygen/bibliography.bib
@@ -22,6 +22,17 @@ @Article{ahlrichs99a
   doi                      = {10.1063/1.480156},
 }
 
+@book{allen2017,
+  title={Computer simulation of liquids},
+  author={Allen, Michael P and Tildesley, Dominic J},
+  year={2017},
+  publisher={Oxford University Press},
+  url={https://global.oup.com/academic/product/computer-simulation-of-liquids-9780198803201},
+  doi={10.1093/oso/9780198803195.001.0001},
+  isbn={9780198803195},
+  edition={2nd},
+}
+
 @Article{andersen83a,
   title                    = {{Rattle: A "velocity" version of the shake algorithm for molecular dynamics calculations}},
   author                   = {Andersen, Hans C.},
@@ -55,6 +66,27 @@ @Article{arnold13c
   issn                     = {1099-4300},
 }
 
+@Article{banchio03a,
+  author  = {Adolfo J. Banchio and John F. Brady},
+  title   = {Accelerated Stokesian dynamics: Brownian motion},
+  journal = {Journal of Chemical Physics},
+  year    = {2003},
+  volume  = {118},
+  number  = {22},
+  pages   = {10323},
+  doi     = {10.1063/1.1571819},
+}
+
+@Article{brady88a,
+  author  = {J. F. Brady and G. Bossis},
+  title   = {Stokesian Dynamics},
+  journal = {Annual Review of Fluid Mechanics},
+  year    = {1988},
+  volume  = {20},
+  pages   = {111-157},
+  doi     = {10.1146/annurev.fl.20.010188.000551},
+}
+
 @Article{brodka04a,
   author    = {Br\'{o}dka, A.},
   title     = {{Ewald summation method with electrostatic layer correction for interactions of point dipoles in slab geometry}},
@@ -91,6 +123,17 @@ @Article{cerda08d
   doi           = {10.1063/1.3000389},
 }
 
+@Article{cortez15a,
+  author  = {Ricardo Cortez and DouglasVarela},
+  title   = {A general system of images for regularized Stokeslets and other elements near a plane wall},
+  journal = {Journal of Computational Physics},
+  year    = {2015},
+  volume  = {285},
+  pages   = {41-54},
+  doi     = {10.1016/j.jcp.2015.01.019},
+  url     = {https://www.sciencedirect.com/science/article/pii/S0021999115000236?via%3Dihub},
+}
+
 @PhdThesis{deserno00b,
   author        = {Deserno, Markus},
   title         = {Counterion condensation for rigid linear polyelectrolytes},
@@ -192,6 +235,17 @@ @Article{dupin08a
 doi = {10.1080/10618560802238242},
 }
 
+@Article{durlofsky87a,
+  author  = {L. Durlofsky and J. F. Brady and G. Bossis},
+  title   = {Dynamic simulation of hydrodynamically interacting particles},
+  journal = {Journal of Fluid Mechanics},
+  year    = {1987},
+  volume  = {180},
+  pages   = {21-49},
+  doi     = {10.1017/S002211208700171X},
+  url     = {https://www.cambridge.org/core/journals/journal-of-fluid-mechanics/article/dynamic-simulation-of-hydrodynamically-interacting-particles/7B6A8C5A04E3D16591FEABE9A8D2778F},
+}
+
 @Article{essmann95a,
   title                    = {{A smooth Particle Mesh Ewald method}},
   author                   = {Essmann, U. and Perera, L. and Berkowitz, M. L. and Darden, T. and Lee, H. and Pedersen, L.},
@@ -265,6 +319,28 @@ @Book{hockney88a
   url       = {https://www.crcpress.com/Computer-Simulation-Using-Particles/Hockney-Eastwood/p/book/9780852743928},
 }
 
+@Article{jancigova16a,
+author = {Jan\v{c}igov\'{a}, Iveta and Cimr\'{a}k, Ivan},
+title = {{Non-uniform force allocation for area preservation in spring network models}},
+journal = {International Journal for Numerical Methods in Biomedical Engineering},
+volume = {32},
+number = {10},
+pages = {e02757},
+year = {2016},
+doi = {10.1002/cnm.2757},
+}
+
+@Article{jeffrey84a,
+  author  = {D. J. Jeffrey and Y. Onishi},
+  title   = {Calculation of the resistance and mobility functions for two unequal rigid spheres in low-Reynolds-number flow},
+  journal = {Journal of Fluid Mechanics},
+  year    = {1984},
+  volume  = {139},
+  pages   = {261-290},
+  doi     = {10.1017/S0022112084000355},
+  url     = {https://www.cambridge.org/core/journals/journal-of-fluid-mechanics/article/calculation-of-the-resistance-and-mobility-functions-for-two-unequal-rigid-spheres-in-lowreynoldsnumber-flow/BDF070924932AD9ABFCDAEC573D1B51A},
+}
+
 @Article{kolafa92a,
   title                    = {{Cutoff errors in the Ewald summation formulae for point charge systems}},
   author                   = {Kolafa, Jiri and Perram, John W.},
@@ -286,17 +362,6 @@ @Book{kruger12a
 doi={10.1007/978-3-8348-2376-2},
 }
 
-@Article{jancigova16a,
-author = {Jan\v{c}igov\'{a}, Iveta and Cimr\'{a}k, Ivan},
-title = {{Non-uniform force allocation for area preservation in spring network models}},
-journal = {International Journal for Numerical Methods in Biomedical Engineering},
-volume = {32},
-number = {10},
-pages = {e02757},
-year = {2016},
-doi = {10.1002/cnm.2757},
-}
-
 @Article{ladd01a,
   title                    = {{Lattice-Boltzmann simulations of particle-fluid suspensions}},
   author                   = {Ladd, A. J. C. and Verberg, R.},
@@ -442,17 +507,6 @@ @Article{sonnenschein85a
 doi = {10.1016/0021-9991(85)90151-2},
 }
 
-@book{allen2017,
-  title={Computer simulation of liquids},
-  author={Allen, Michael P and Tildesley, Dominic J},
-  year={2017},
-  publisher={Oxford University Press},
-  url={https://global.oup.com/academic/product/computer-simulation-of-liquids-9780198803201},
-  doi={10.1093/oso/9780198803195.001.0001},
-  isbn={9780198803195},
-  edition={2nd},
-}
-
 @Article{swope92a,
 author = {Swope, William C. and Ferguson, David M.},
 title = {{Alternative expressions for energies and forces due to angle bending and torsional energy}},

diff --git a/doc/sphinx/CMakeLists.txt b/doc/sphinx/CMakeLists.txt
@@ -43,6 +43,7 @@ if(SPHINX_FOUND)
       "${CMAKE_CURRENT_SOURCE_DIR}/magnetostatics.rst"
       "${CMAKE_CURRENT_SOURCE_DIR}/particles.rst"
       "${CMAKE_CURRENT_SOURCE_DIR}/running.rst"
+      "${CMAKE_CURRENT_SOURCE_DIR}/sd.rst"
       "${CMAKE_CURRENT_SOURCE_DIR}/system_manipulation.rst"
       "${CMAKE_CURRENT_SOURCE_DIR}/system_setup.rst"
       "${CMAKE_CURRENT_SOURCE_DIR}/ug.rst"

diff --git a/doc/sphinx/conf.py.in b/doc/sphinx/conf.py.in
@@ -256,7 +256,12 @@ if sphinx.version_info[:3] < (2, 0, 1):
     autodoc_mock_imports.append('MDAnalysis')
 
 def setup(app):
-    app.add_css_file('bibtex.css')
+
+    # Function name was changed from before version 1.8
+    try:
+        app.add_css_file('bibtex.css')
+    except AttributeError:
+        app.add_stylesheet('bibtex.css')
 
 
 # -- Options for Cython ---------------------------------------------------

diff --git a/doc/sphinx/sd.rst b/doc/sphinx/sd.rst
@@ -0,0 +1,69 @@
+.. _Stokesian_Dynamics:
+
+Stokesian Dynamics
+==================
+
+The Stokesian Dynamics method allows to study the behaviour of spherical
+particles in a viscous fluid. It is targeted at systems with very low Reynolds
+numbers. In such systems, particles stop moving immediately as soon as any
+force on them is removed. In other words, motion has no memory of the past. 
+
+The integration scheme is relatively simple. Only the particle's positions,
+radii and forces (including torques) are needed to compute the momentary
+velocities (including angular velocities). The particle positions are
+integrated by simple Euler scheme.
+
+The Stokesian Dynamics method is outlined in :cite:`durlofsky87a`
+
+..note:: The computation of the velocities is an approximation with good results in the far field.
+
+..note:: The Stokesian Dynamics method is only available for open systems, i.e. no periodic boundary conditions are supported. The box size has no effect either.
+
+.. _Setting up an SD system:
+
+Setting up an SD system:
+------------------------
+
+The following minimal example illustrates how to use the SDM in |es|::
+
+    import espressomd
+    system = espressomd.System(box_l=[1.0, 1.0, 1.0])
+    system.time_step = 0.01
+    system.cell_system.skin = 0.4
+    system.part.add(pos=[0, 0, 0], rotation=[1, 0, 0])
+    system.integrator.set_sd(viscosity=1.0, radii={0: 1.0})
+
+    system.integrator.run(100)
+
+Because there is no force on the particle yet, nothing will move. You will need
+to add your own actors to the system. The parameter ``radii`` is a dictionary
+that maps particle types to different radii. ``viscosity`` is the dynamic
+viscosity of the ambient infinite fluid. There are additional optional
+parameters for ``set_sd()``. For more information, see
+:py:meth:`espressomd.integrate.IntegratorHandle.set_sd()`.
+
+Note that this setup represents a system at zero temperature. In order to
+thermalize the system, the SD thermostat can be activated via an additional
+command::
+
+    system.thermostat.set_sd(kT=1.0, seed=43)
+
+``kT`` denotes the desired temperature of the system, and ``seed`` the seed for
+the random number generator of the Stokesian Dynamics thermostat. It is
+independent from other random number generators inside |es|.
+
+.. _Important_SD:
+
+Important:
+----------
+
+The particles must be prevented from overlapping. It is mathematically allowed for the particles to overlap to a certain degree. However, once the distance of the sphere centers is less than 2/3 of the sphere diameter, the mobility matrix is no longer positive definite and the Stokesian Dynamics integrator will fail. Therefore, the particle centers must be kept apart from each other by a strongly repulsive potential, for example the WCA potential that is set to the appropriate particle radius (for more information about the available interaction types see :ref:`Non-bonded interactions`).
+
+The current implementation of SD only includes the far field approximation.
+The near field (so-called lubrication) correction is planned. For now,
+Stokesian Dynamics provides a good approximation of the hydrodynamics
+in dilute systems where the average distance between particles is several
+sphere diameters.
+
+
+.. check math claims
diff --git a/doc/sphinx/ug.rst b/doc/sphinx/ug.rst
@@ -14,6 +14,7 @@
     system_manipulation.rst
     constraints.rst
     lb.rst
+    sd.rst
     analysis.rst
     io.rst
     visualization.rst

diff --git a/doc/sphinx/zrefs.bib b/doc/sphinx/zrefs.bib
@@ -1,7 +1,6 @@
-% This file was created with JabRef 2.9.
 % Encoding: US-ASCII
 @article{ahlrichs99,
-author = {Ahlrichs,Patrick  and Dünweg,Burkhard },
+author = {Ahlrichs,Patrick  and D\"{u}nweg,Burkhard },
 title = {Simulation of a single polymer chain in solution by combining lattice Boltzmann and molecular dynamics},
 journal = {The Journal of Chemical Physics},
 volume = {111},
@@ -1075,3 +1074,16 @@ @article{yaghoubi2015a
   year={2015},
   publisher={IOP Publishing}
 }
+
+@Article{durlofsky87a,
+  author  = {L. Durlofsky and J. F. Brady and G. Bossis},
+  title   = {Dynamic simulation of hydrodynamically interacting particles},
+  journal = {Journal of Fluid Mechanics},
+  year    = {1987},
+  volume  = {180},
+  pages   = {21-49},
+  doi     = {10.1017/S002211208700171X},
+  url     = {https://www.cambridge.org/core/journals/journal-of-fluid-mechanics/article/dynamic-simulation-of-hydrodynamically-interacting-particles/7B6A8C5A04E3D16591FEABE9A8D2778F},
+}
+
+@Comment{jabref-meta: databaseType:bibtex;}
diff --git a/libs/CMakeLists.txt b/libs/CMakeLists.txt
@@ -4,3 +4,7 @@ target_compile_definitions(Random123 INTERFACE R123_USE_MULHILO64_C99)
 
 add_library(h5xx INTERFACE)
 target_include_directories(h5xx SYSTEM INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/h5xx)
+
+if(STOKESIAN_DYNAMICS OR STOKESIAN_DYNAMICS_GPU)
+  add_subdirectory(stokesian_dynamics)
+endif()
diff --git a/libs/stokesian_dynamics/CMakeLists.txt b/libs/stokesian_dynamics/CMakeLists.txt
@@ -0,0 +1,36 @@
+if(STOKESIAN_DYNAMICS_GPU)
+  add_gpu_library(sd_gpu SHARED sd_gpu.cu)
+  target_link_libraries(sd_gpu PRIVATE
+    ${CUDA_CUBLAS_LIBRARIES}
+    ${CUDA_cusolver_LIBRARY})
+  target_compile_definitions(sd_gpu PRIVATE SD_USE_THRUST SD_USE_GPU)
+  if(HIP)
+    target_compile_definitions(sd_gpu PRIVATE SD_USE_ROCM)
+  else()
+    target_compile_definitions(sd_gpu PRIVATE SD_USE_CUDA)
+  endif()
+  install(TARGETS sd_gpu LIBRARY DESTINATION ${PYTHON_INSTDIR}/espressomd)
+endif()
+
+if(STOKESIAN_DYNAMICS)
+  add_library(sd_cpu SHARED sd_cpu.cpp)
+  target_link_libraries(sd_cpu 
+    PRIVATE
+      ${BLAS_LIBRARIES}
+      ${LAPACK_LIBRARIES}
+      Boost::boost
+      Random123)
+
+  # In case the GPU is used, Thrust is present and can be used to parallelize
+  # the CPU code, too. The standard compiler needs to be told the location of
+  # Thrust
+  if(STOKESIAN_DYNAMICS_GPU)
+    target_compile_definitions(sd_cpu PUBLIC SD_USE_THRUST PRIVATE THRUST_DEVICE_SYSTEM=4)
+    #target_compile_definitions(sd_cpu PUBLIC USE_THRUST)
+    target_include_directories(sd_cpu PRIVATE 
+      ${CUDA_INCLUDE_DIRS}
+      ${HIP_ROOT_DIR}/include
+      ${ROCM_HOME}/include)
+  endif()
+  install(TARGETS sd_cpu LIBRARY DESTINATION ${PYTHON_INSTDIR}/espressomd)
+endif()