diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 6fc32ce34e2..5cfe16c18e9 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -77,7 +77,8 @@ default: CC: 'gcc-9' CXX: 'g++-9' script: - - export with_cuda=false myconfig=default with_coverage=true with_scafacos=true + - export with_cuda=false myconfig=default with_coverage=true + - export with_scafacos=true with_stokesian_dynamics=true - bash maintainer/CI/build_cmake.sh tags: - docker @@ -90,7 +91,8 @@ maxset: CC: 'gcc-9' CXX: 'g++-9' script: - - export with_cuda=false myconfig=maxset with_coverage=true with_scafacos=true + - export with_cuda=false myconfig=maxset with_coverage=true + - export with_scafacos=true with_stokesian_dynamics=true - bash maintainer/CI/build_cmake.sh tags: - docker @@ -125,20 +127,20 @@ ubuntu:wo-dependencies: debian:10: <<: *global_job_definition stage: build - image: docker.pkg.github.com/espressomd/docker/debian:446ff604bbfa63f30ddb462697fa0d0dc2630460 + image: docker.pkg.github.com/espressomd/docker/debian:d496478230db4e5c286680e3bdc1621af1fccffc script: - - export with_cuda=false myconfig=maxset make_check_python=false + - export with_cuda=false myconfig=maxset make_check_python=false with_stokesian_dynamics=true - bash maintainer/CI/build_cmake.sh tags: - docker - linux -opensuse:15.1: +opensuse:15.2: <<: *global_job_definition stage: build - image: docker.pkg.github.com/espressomd/docker/opensuse:446ff604bbfa63f30ddb462697fa0d0dc2630460 + image: docker.pkg.github.com/espressomd/docker/opensuse:d496478230db4e5c286680e3bdc1621af1fccffc script: - - export with_cuda=false myconfig=maxset make_check_python=false + - export with_cuda=false myconfig=maxset make_check_python=false with_stokesian_dynamics=true - bash maintainer/CI/build_cmake.sh tags: - docker @@ -149,7 +151,7 @@ centos:7: stage: build image: docker.pkg.github.com/espressomd/docker/centos:446ff604bbfa63f30ddb462697fa0d0dc2630460 script: - - export with_cuda=false myconfig=maxset make_check_python=true + - export with_cuda=false myconfig=maxset make_check_python=true with_stokesian_dynamics=true - bash maintainer/CI/build_cmake.sh tags: - docker @@ -160,7 +162,7 @@ fedora:32: stage: build image: docker.pkg.github.com/espressomd/docker/fedora:fc7628d32de0fce605976ba9edebe7eff186e618 script: - - export with_cuda=false myconfig=maxset make_check_python=false + - export with_cuda=false myconfig=maxset make_check_python=false with_stokesian_dynamics=true - bash maintainer/CI/build_cmake.sh tags: - docker @@ -187,8 +189,9 @@ clang-sanitizer: CC: 'clang-9' CXX: 'clang++-9' script: - - export myconfig=maxset with_cuda=true with_cuda_compiler=clang with_coverage=false with_scafacos=true + - export myconfig=maxset with_cuda=true with_cuda_compiler=clang with_coverage=false - export with_static_analysis=true test_timeout=900 with_asan=true with_ubsan=true + - export with_scafacos=true with_stokesian_dynamics=true - bash maintainer/CI/build_cmake.sh timeout: 2h tags: @@ -203,7 +206,8 @@ cuda10-maxset: CC: 'gcc-8' CXX: 'g++-8' script: - - export myconfig=maxset with_cuda=true with_coverage=false with_scafacos=true test_timeout=900 srcdir=${CI_PROJECT_DIR} + - export myconfig=maxset with_cuda=true with_coverage=false test_timeout=900 srcdir=${CI_PROJECT_DIR} + - export with_scafacos=true with_stokesian_dynamics=true - bash maintainer/CI/build_cmake.sh artifacts: paths: @@ -217,13 +221,14 @@ cuda10-maxset: cuda9-maxset: <<: *global_job_definition stage: build - image: docker.pkg.github.com/espressomd/docker/ubuntu-18.04:06b6216c7aa3555bcf28c90734dbb84e7285c96f + image: docker.pkg.github.com/espressomd/docker/ubuntu-18.04:d496478230db4e5c286680e3bdc1621af1fccffc variables: CC: 'gcc-6' CXX: 'g++-6' GCOV: 'gcov-6' script: - - export myconfig=maxset with_cuda=true with_coverage=true with_scafacos=true test_timeout=900 srcdir=${CI_PROJECT_DIR} + - export myconfig=maxset with_cuda=true with_coverage=true test_timeout=900 srcdir=${CI_PROJECT_DIR} + - export with_scafacos=true with_stokesian_dynamics=true - bash maintainer/CI/build_cmake.sh artifacts: paths: @@ -307,7 +312,8 @@ installation: CC: 'gcc-8' CXX: 'g++-8' script: - - export myconfig=maxset with_cuda=true with_coverage=false with_scafacos=true make_check_unit_tests=false make_check_python=false + - export myconfig=maxset with_cuda=true with_coverage=false make_check_unit_tests=false make_check_python=false + - export with_scafacos=true with_stokesian_dynamics=true - export srcdir=${CI_PROJECT_DIR} build_type=Release - bash maintainer/CI/build_cmake.sh - cd build @@ -335,7 +341,8 @@ empty: CC: 'clang-9' CXX: 'clang++-9' script: - - export myconfig=empty with_cuda=true with_cuda_compiler=clang with_static_analysis=true with_scafacos=false + - export myconfig=empty with_cuda=true with_cuda_compiler=clang with_static_analysis=true + - export with_scafacos=false with_stokesian_dynamics=false - bash maintainer/CI/build_cmake.sh tags: - docker @@ -347,9 +354,10 @@ empty: rocm-maxset: <<: *global_job_definition stage: build - image: docker.pkg.github.com/espressomd/docker/rocm:06b6216c7aa3555bcf28c90734dbb84e7285c96f + image: docker.pkg.github.com/espressomd/docker/rocm:d496478230db4e5c286680e3bdc1621af1fccffc script: - export myconfig=maxset with_cuda=true with_cuda_compiler=hip + - export with_stokesian_dynamics=false - bash maintainer/CI/build_cmake.sh tags: - amdgpu @@ -360,6 +368,7 @@ rocm:latest: image: docker.pkg.github.com/espressomd/docker/rocm:latest_base script: - export myconfig=maxset with_cuda=true with_cuda_compiler=hip + - export with_stokesian_dynamics=false - bash maintainer/CI/build_cmake.sh tags: - amdgpu diff --git a/AUTHORS b/AUTHORS index 0bfa201ff11..72f4763c3be 100644 --- a/AUTHORS +++ b/AUTHORS @@ -39,6 +39,7 @@ Ashreya Jayaram Ben Reynwar (formerly: Reynolds) Bogdan Tanygin Cameron Stewart +Carl Georg Biermann Christian Haege Christoph Lohrmann Christoph Schneider diff --git a/CMakeLists.txt b/CMakeLists.txt index fde9eaf1344..37f6a1743a5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -74,6 +74,7 @@ option(WITH_CUDA "Build with GPU support" OFF) option_if_available(WITH_HDF5 "Build with HDF5 support" ON) option(WITH_TESTS "Enable tests" ON) option_if_available(WITH_SCAFACOS "Build with ScaFaCoS support" OFF) +option_if_available(WITH_STOKESIAN_DYNAMICS "Build with Stokesian Dynamics" ON) option(WITH_BENCHMARKS "Enable benchmarks" OFF) option(WITH_VALGRIND_INSTRUMENTATION "Build with valgrind instrumentation markers" OFF) @@ -253,6 +254,62 @@ 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(WITH_STOKESIAN_DYNAMICS) + if(BLAS AND LAPACK) + set(STOKESIAN_DYNAMICS 1) + endif() + if(CUDA) + set(STOKESIAN_DYNAMICS_GPU 1) + if(HIP) + find_gpu_library(VARNAME CUDA_CUBLAS_LIBRARIES NAMES rocblas) + find_gpu_library(VARNAME CUDA_cusolver_LIBRARY NAMES rocsolver) + else() + if(CMAKE_CUDA_COMPILER_VERSION VERSION_GREATER_EQUAL 10.0) + # when importing package cublas, CMake versions older than 3.12.2 will + # look for a library cublas_device that no longer exists in CUDA 10 + cmake_minimum_required(VERSION 3.12.2) + endif() + find_gpu_library(VARNAME CUDA_CUBLAS_LIBRARIES NAMES cublas) + find_gpu_library(VARNAME CUDA_cusolver_LIBRARY NAMES cusolver) + endif() + endif(CUDA) + if(NOT (STOKESIAN_DYNAMICS + OR (STOKESIAN_DYNAMICS_GPU AND CUDA_CUBLAS_LIBRARIES + AND CUDA_cusolver_LIBRARY)) + AND NOT WITH_STOKESIAN_DYNAMICS_IS_DEFAULT_VALUE) + message( + FATAL_ERROR + "Optional feature Stokesian Dynamics explicitly requested, but dependencies not found." + ) + endif() +endif(WITH_STOKESIAN_DYNAMICS) + +if(WITH_STOKESIAN_DYNAMICS) + set(CMAKE_INSTALL_LIBDIR + "${CMAKE_INSTALL_PREFIX}/${PYTHON_INSTDIR}/espressomd") + cmake_minimum_required(VERSION 3.11) + include(FetchContent) + FetchContent_Declare( + stokesian_dynamics + GIT_REPOSITORY https://github.com/hmenke/espresso-stokesian-dynamics.git + GIT_TAG 860ceac208175) + FetchContent_GetProperties(stokesian_dynamics) + if(NOT stokesian_dynamics_POPULATED) + FetchContent_Populate(stokesian_dynamics) + add_subdirectory(${stokesian_dynamics_SOURCE_DIR} + ${stokesian_dynamics_BINARY_DIR}) + endif() +endif(WITH_STOKESIAN_DYNAMICS) + if(WITH_VALGRIND_INSTRUMENTATION) find_package(PkgConfig) pkg_check_modules(VALGRIND valgrind) diff --git a/cmake/cmake_config.cmakein b/cmake/cmake_config.cmakein index 9389e251138..1bc7dcd59e2 100644 --- a/cmake/cmake_config.cmakein +++ b/cmake/cmake_config.cmakein @@ -13,6 +13,14 @@ #cmakedefine GSL +#cmakedefine BLAS + +#cmakedefine LAPACK + +#cmakedefine STOKESIAN_DYNAMICS + +#cmakedefine STOKESIAN_DYNAMICS_GPU + #cmakedefine VALGRIND_INSTRUMENTATION #define PACKAGE_NAME "${PROJECT_NAME}" diff --git a/doc/doxygen/bibliography.bib b/doc/doxygen/bibliography.bib index 967580da4d2..6df6aa756c6 100644 --- a/doc/doxygen/bibliography.bib +++ b/doc/doxygen/bibliography.bib @@ -22,6 +22,16 @@ @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}, + 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 +65,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 +122,16 @@ @Article{cerda08d doi = {10.1063/1.3000389}, } +@Article{cortez15a, + author = {Cortez, Ricardo and Varela, Douglas}, + 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}, +} + @PhdThesis{deserno00b, author = {Deserno, Markus}, title = {Counterion condensation for rigid linear polyelectrolytes}, @@ -192,6 +233,16 @@ @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}, +} + @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 +316,27 @@ @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}, +} + @Article{kolafa92a, title = {{Cutoff errors in the Ewald summation formulae for point charge systems}}, author = {Kolafa, Jiri and Perram, John W.}, @@ -286,17 +358,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.}, @@ -329,6 +390,17 @@ @Article{marsili10a year = {2010}, } +@article{martys99a, + title = {{Velocity Verlet algorithm for dissipative-particle-dynamics-based models of suspensions}}, + author = {Martys, Nicos S. and Mountain, Raymond D.}, + journal = {Physical Review E}, + volume = {59}, + number = {3}, + pages = {3733--3736}, + year = {1999}, + doi = {10.1103/PhysRevE.59.3733}, +} + @Article{moss96a, author={Moss, G. P.}, journal={Pure and Applied Chemistry}, @@ -351,6 +423,17 @@ @Article{neumann85b doi = {10.1063/1.448553}, } +@article{omelyan98a, +author = {Omelyan, Igor P.}, +title = {{On the numerical integration of motion for rigid polyatomics: The modified quaternion approach}}, +journal = {Computers in Physics}, +volume = {12}, +number = {1}, +pages = {97-103}, +year = {1998}, +doi = {10.1063/1.168642}, +} + @article{panneton06a, author = {Panneton, Fran\c{c}ois and L'Ecuyer, Pierre and Matsumoto, Makoto}, title = {Improved long-period generators based on linear recurrences modulo 2}, @@ -442,38 +525,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}, - doi={10.1093/oso/9780198803195.001.0001}, - isbn={9780198803195}, - edition={2nd}, -} - -@article{martys99a, - title = {{Velocity Verlet algorithm for dissipative-particle-dynamics-based models of suspensions}}, - author = {Martys, Nicos S. and Mountain, Raymond D.}, - journal = {Physical Review E}, - volume = {59}, - number = {3}, - pages = {3733--3736}, - year = {1999}, - doi = {10.1103/PhysRevE.59.3733}, -} - -@article{omelyan98a, -author = {Omelyan, Igor P.}, -title = {{On the numerical integration of motion for rigid polyatomics: The modified quaternion approach}}, -journal = {Computers in Physics}, -volume = {12}, -number = {1}, -pages = {97-103}, -year = {1998}, -doi = {10.1063/1.168642}, -} - @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/installation.rst b/doc/sphinx/installation.rst index 0644d525d8d..324d0214739 100644 --- a/doc/sphinx/installation.rst +++ b/doc/sphinx/installation.rst @@ -485,9 +485,15 @@ Finally, there is a flag for debugging: Features marked as experimental ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -Some of the above features are marked as EXPERIMENTAL. Activating these features can have unexpected side effects and some of them have known issues. If you activate any of these features, you should understand the corresponding source code and do extensive testing. Furthermore, it is necessary to define ``EXPERIMENTAL_FEATURES`` in :file:`myconfig.hpp`. +Some of the above features are marked as EXPERIMENTAL. Activating these +features can have unexpected side effects and some of them have known issues. +If you activate any of these features, you should understand the corresponding +source code and do extensive testing. Furthermore, it is necessary to define +``EXPERIMENTAL_FEATURES`` in :file:`myconfig.hpp`. +.. _External features: + External features ^^^^^^^^^^^^^^^^^ @@ -508,6 +514,14 @@ using a CMake flag (see :ref:`Options and Variables`). - ``GSL`` Enables features relying on the GNU Scientific Library, e.g. :meth:`espressomd.cluster_analysis.Cluster.fractal_dimension`. +- ``STOKESIAN_DYNAMICS`` Enables the Stokesian Dynamics feature for CPU + (see :ref:`Stokesian Dynamics`). Requires BLAS and LAPACK. + +- ``STOKESIAN_DYNAMICS_GPU`` Enables the Stokesian Dynamics feature for GPU + (see :ref:`Stokesian Dynamics`). Requires thrust/cuBLAS/cuSolver for NVIDIA + GPUs or rocrand/rocthrust/rocblas/rocsolver for AMD GPUs. + Requires ``EXPERIMENTAL_FEATURES``. + .. _Configuring: @@ -563,11 +577,11 @@ Then you can simply compile two different versions of |es| via: .. code-block:: bash - cd builddir1 + cd $builddir1 cmake .. make - cd builddir2 + cd $builddir2 cmake .. make @@ -670,6 +684,8 @@ options are available: * ``WITH_SCAFACOS``: Build with ScaFaCoS support +* ``WITH_STOKESIAN_DYNAMICS`` Build with Stokesian Dynamics support + * ``WITH_VALGRIND_INSTRUMENTATION``: Build with valgrind instrumentation markers @@ -697,6 +713,40 @@ correct one with ``CUDA_BIN_PATH=/usr/local/cuda-10.0 cmake .. -DWITH_CUDA=ON`` path with ``-DCMAKE_CXX_FLAGS=--cuda-path=/usr/local/cuda-10.0``). +.. _Configuring without a network connection: + +Configuring without a network connection +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Several :ref:`external features ` in |es| rely on +external libraries that are downloaded automatically by CMake. When a +network connection cannot be established due to firewall restrictions, +the CMake logic needs editing: + +* ``WITH_HDF5``: when cloning |es|, the :file:`/libs/h5xx` folder will be + a git submodule containing a :file:`.git` subfolder. To prevent CMake from + updating this submodule with git, delete the corresponding command with: + + .. code-block:: bash + + sed -i '/execute_process(COMMAND ${GIT_EXECUTABLE} submodule update -- libs\/h5xx/,+1 d' CMakeLists.txt + + When installing a release version of |es|, no network communication + is needed for HDF5. + +* ``WITH_STOKESIAN_DYNAMICS``: this library is installed using `FetchContent + `_. + The repository URL can be found in the ``GIT_REPOSITORY`` field of the + corresponding ``FetchContent_Declare()`` command. The ``GIT_TAG`` field + provides the commit. Clone this repository locally next to the |es| + folder and edit the |es| build system such that ``GIT_REPOSITORY`` points + to the absolute path of the Stokesian Dynamics clone, for example with: + + .. code-block:: bash + + sed -ri 's|GIT_REPOSITORY +.+stokesian-dynamics.git|GIT_REPOSITORY /work/username/stokesian_dynamics|' CMakeLists.txt + + Compiling, testing and installing --------------------------------- diff --git a/doc/sphinx/introduction.rst b/doc/sphinx/introduction.rst index 79b1adbe869..552cbaf72e5 100644 --- a/doc/sphinx/introduction.rst +++ b/doc/sphinx/introduction.rst @@ -435,6 +435,10 @@ report so to the developers using the instructions in :ref:`Contributing`. +--------------------------------+------------------------+------------------+------------+ | Quaternion Integrator | Core | Good | Yes | +--------------------------------+------------------------+------------------+------------+ +| Stokesian Dynamics on CPU | Single | None | Yes | ++--------------------------------+------------------------+------------------+------------+ +| Stokesian Dynamics on GPU | Experimental | None | No | ++--------------------------------+------------------------+------------------+------------+ | **Interactions** | +--------------------------------+------------------------+------------------+------------+ | Short-range Interactions | Core | Core | Yes | diff --git a/doc/sphinx/running.rst b/doc/sphinx/running.rst index ea6b9367ce2..aa2930dbe3f 100644 --- a/doc/sphinx/running.rst +++ b/doc/sphinx/running.rst @@ -273,3 +273,75 @@ This convenience function only exists for historical reasons and remains availab for backward compatibility. New scripts should setup the steepest descent integrator with the :meth:`~espressomd.integrate.IntegratorHandle.set_steepest_descent` handle directly. + +.. _Stokesian Dynamics: + +Stokesian Dynamics +------------------ + +.. note:: + + Requires ``STOKESIAN_DYNAMICS`` external feature, enabled with + ``-DWITH_STOKESIAN_DYNAMICS=ON``. + +:meth:`espressomd.integrate.IntegratorHandle.set_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 come to a rest almost 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 particles' positions, +radii and forces (including torques) are needed to compute the momentary +velocities (including angular velocities). The particle positions are +integrated by the simple Euler scheme. + +The computation of the velocities is an approximation with good results +in the far field. +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. + +The Stokesian Dynamics method is outlined in :cite:`durlofsky87a`. + +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_stokesian_dynamics(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_stokesian_dynamics()``. For more information, see +:py:meth:`espressomd.integrate.IntegratorHandle.set_stokesian_dynamics()`. + +Note that this setup represents a system at zero temperature. In order to +thermalize the system, the SD thermostat needs to be activated (see +:ref:`Stokesian thermostat`). + +.. _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. diff --git a/doc/sphinx/system_setup.rst b/doc/sphinx/system_setup.rst index be5d738e122..bafbafdd0df 100644 --- a/doc/sphinx/system_setup.rst +++ b/doc/sphinx/system_setup.rst @@ -473,6 +473,26 @@ can be omitted in subsequent calls of ``set_brownian()``. It is the user's responsibility to decide whether the thermostat should be deterministic (by using a fixed seed) or not (by using a randomized seed). +.. _Stokesian thermostat: + +Stokesian thermostat +~~~~~~~~~~~~~~~~~~~~ + +.. note:: + + Requires ``STOKESIAN_DYNAMICS`` external feature, enabled with + ``-DWITH_STOKESIAN_DYNAMICS=ON``. + +In order to thermalize a Stokesian Dynamics simulation, the SD thermostat +needs to be activated via:: + + system.thermostat.set_stokesian(kT=1.0, seed=43) + +where ``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 the other random number generators in |es|. + + .. _CUDA: CUDA diff --git a/doc/sphinx/zrefs.bib b/doc/sphinx/zrefs.bib index 132aab7b3e1..e085d956367 100644 --- 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 {B}oltzmann and molecular dynamics}, journal = {The Journal of Chemical Physics}, volume = {111}, @@ -998,3 +997,13 @@ @article{yaghoubi2015a doi={10.1209/0295-5075/110/24002}, 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}, +} diff --git a/libs/CMakeLists.txt b/libs/CMakeLists.txt index ac344dca2fb..2e876e3dbb3 100644 --- a/libs/CMakeLists.txt +++ b/libs/CMakeLists.txt @@ -1,6 +1,8 @@ add_library(Random123 INTERFACE) -target_include_directories(Random123 SYSTEM INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/Random123-1.09/include) +target_include_directories( + Random123 SYSTEM INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/Random123-1.09/include) target_compile_definitions(Random123 INTERFACE R123_USE_MULHILO64_C99) add_library(h5xx INTERFACE) -target_include_directories(h5xx SYSTEM INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/h5xx) +target_include_directories(h5xx SYSTEM + INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/h5xx) diff --git a/maintainer/CI/build_cmake.sh b/maintainer/CI/build_cmake.sh index 341d87f63ae..fa351b7ac82 100755 --- a/maintainer/CI/build_cmake.sh +++ b/maintainer/CI/build_cmake.sh @@ -97,6 +97,7 @@ set_default_value with_cuda_compiler "nvcc" set_default_value build_type "RelWithAssert" set_default_value with_ccache false set_default_value with_scafacos false +set_default_value with_stokesian_dynamics false set_default_value test_timeout 300 set_default_value hide_gpu false @@ -120,6 +121,13 @@ if [ "${with_ccache}" = true ]; then fi if [ "${with_scafacos}" = true ]; then cmake_params="${cmake_params} -DWITH_SCAFACOS=ON" +else + cmake_params="${cmake_params} -DWITH_SCAFACOS=OFF" +fi +if [ "${with_stokesian_dynamics}" = true ]; then + cmake_params="${cmake_params} -DWITH_STOKESIAN_DYNAMICS=ON" +else + cmake_params="${cmake_params} -DWITH_STOKESIAN_DYNAMICS=OFF" fi if [ "${with_fftw}" = true ]; then diff --git a/pylint.log b/pylint.log new file mode 100644 index 00000000000..e69de29bb2d diff --git a/samples/dancing.py b/samples/dancing.py new file mode 100644 index 00000000000..67c0b5390b9 --- /dev/null +++ b/samples/dancing.py @@ -0,0 +1,101 @@ +# +# Copyright (C) 2019-2020 The ESPResSo project +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +""" +Stokesian Dynamics simulation of particle sedimentation. +Reproduce the trajectory in Figure 5b from :cite:`durlofsky87a`. +""" +import espressomd +import espressomd.constraints +import espressomd.observables +import espressomd.accumulators +import numpy as np +import matplotlib.pyplot as plt + +import argparse + +parser = argparse.ArgumentParser(epilog=__doc__) +group = parser.add_mutually_exclusive_group() +group.add_argument('--cpu', action='store_true', help='Use CPU implementation') +group.add_argument('--gpu', action='store_true', help='Use GPU implementation') +group = parser.add_mutually_exclusive_group() +group.add_argument('--ft', action='store_true', help='Use FT approximation') +group.add_argument('--fts', action='store_true', help='Use FTS approximation') +args = parser.parse_args() + +required_features = ["STOKESIAN_DYNAMICS"] + +if args.gpu: + print("Using GPU implementation") + required_features.append("CUDA") + required_features.append("STOKESIAN_DYNAMICS_GPU") + sd_device = "gpu" +else: + print("Using CPU implementation") + sd_device = "cpu" + +if args.ft: + print("Using FT approximation method") + sd_method = "ft" +else: + print("Using FTS approximation method") + sd_method = "fts" + +espressomd.assert_features(required_features) + +system = espressomd.System(box_l=[10, 10, 10]) +system.time_step = 1.5 +system.cell_system.skin = 0.4 +system.periodicity = [False, False, False] + +system.integrator.set_stokesian_dynamics( + viscosity=1.0, radii={0: 1.0}, device=sd_device, + approximation_method=sd_method) + +system.part.add(pos=[-5, 0, 0], rotation=[1, 1, 1]) +system.part.add(pos=[0, 0, 0], rotation=[1, 1, 1]) +system.part.add(pos=[7, 0, 0], rotation=[1, 1, 1]) + +gravity = espressomd.constraints.Gravity(g=[0, -1, 0]) +system.constraints.add(gravity) + +obs = espressomd.observables.ParticlePositions(ids=system.part[:].id) +acc = espressomd.accumulators.TimeSeries(obs=obs, delta_N=1) +system.auto_update_accumulators.add(acc) +acc.update() +intsteps = int(10500 / system.time_step) +system.integrator.run(intsteps) + +positions = acc.time_series() +ref_data = "../testsuite/python/data/dancing.txt" +data = np.loadtxt(ref_data) + +for i in range(3): + plt.plot(positions[:, i, 0], positions[:, i, 1], linestyle='solid') + +plt.gca().set_prop_cycle(None) + +for i in range(0, 6, 2): + plt.plot(data[:, i], data[:, i + 1], linestyle='dashed') + +plt.title("Trajectory of sedimenting spheres\nsolid line: simulation " + "({} on {}), dashed line: paper (FTS)" + .format(sd_method.upper(), sd_device.upper())) +plt.xlabel("x") +plt.ylabel("y") +plt.show() diff --git a/src/config/features.def b/src/config/features.def index f72b4f287f0..a18167fafa2 100644 --- a/src/config/features.def +++ b/src/config/features.def @@ -24,6 +24,12 @@ COLLISION_DETECTION NPT ENGINE implies ROTATION, EXTERNAL_FORCES PARTICLE_ANISOTROPY implies ROTATION +STOKESIAN_DYNAMICS requires BLAS and LAPACK +STOKESIAN_DYNAMICS implies ROTATION +# Needs cuBLAS, which is expected to be included with CUDA +STOKESIAN_DYNAMICS_GPU requires CUDA +STOKESIAN_DYNAMICS_GPU implies ROTATION +STOKESIAN_DYNAMICS_GPU requires EXPERIMENTAL_FEATURES /* Rotation */ ROTATION @@ -96,8 +102,13 @@ ADDITIONAL_CHECKS # External switches # Switches that are set by configure or gcc, not to be set manually +# All these switches must also be present in cmake/cmake_config.cmakein CUDA external FFTW external H5MD external SCAFACOS external GSL external +BLAS external +LAPACK external +STOKESIAN_DYNAMICS external +STOKESIAN_DYNAMICS_GPU external diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 38c7e743810..034e3f4468f 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -109,3 +109,13 @@ add_subdirectory(virtual_sites) if(WITH_UNIT_TESTS) add_subdirectory(unit_tests) endif(WITH_UNIT_TESTS) + +if(STOKESIAN_DYNAMICS OR STOKESIAN_DYNAMICS_GPU) + add_subdirectory(stokesian_dynamics) + if(STOKESIAN_DYNAMICS) + target_link_libraries(EspressoCore PRIVATE StokesianDynamics::sd_cpu) + endif() + if(STOKESIAN_DYNAMICS_GPU) + target_link_libraries(EspressoCore PRIVATE StokesianDynamics::sd_gpu) + endif() +endif() diff --git a/src/core/communication.cpp b/src/core/communication.cpp index 7e0be01d886..8dc7a8b07e5 100644 --- a/src/core/communication.cpp +++ b/src/core/communication.cpp @@ -57,6 +57,7 @@ #include "particle_data.hpp" #include "pressure.hpp" #include "rotation.hpp" +#include "stokesian_dynamics/sd_interface.hpp" #include "virtual_sites.hpp" #include "electrostatics_magnetostatics/coulomb.hpp" diff --git a/src/core/event.cpp b/src/core/event.cpp index e89f7410001..15395fd6418 100644 --- a/src/core/event.cpp +++ b/src/core/event.cpp @@ -287,6 +287,15 @@ void on_parameter_change(int field) { Scafacos::update_system_params(); } #endif +#endif +#ifdef STOKESIAN_DYNAMICS + if (integ_switch == INTEG_METHOD_SD) { + if (box_geo.periodic(0) || box_geo.periodic(1) || box_geo.periodic(2)) + runtimeErrorMsg() << "Illegal box periodicity for Stokesian Dynamics: " + << box_geo.periodic(0) << " " << box_geo.periodic(1) + << " " << box_geo.periodic(2) << "\n" + << " Required: 0 0 0\n"; + } #endif case FIELD_MIN_GLOBAL_CUT: case FIELD_SKIN: { diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index 1687264a98d..c261d477943 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -40,6 +40,7 @@ #include "event.hpp" #include "forces.hpp" #include "global.hpp" +#include "grid.hpp" #include "grid_based_algorithms/lb_interface.hpp" #include "grid_based_algorithms/lb_particle_coupling.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" @@ -53,6 +54,7 @@ #include "integrators/brownian_inline.hpp" #include "integrators/langevin_inline.hpp" #include "integrators/steepest_descent.hpp" +#include "integrators/stokesian_dynamics_inline.hpp" #include "integrators/velocity_verlet_inline.hpp" #include "integrators/velocity_verlet_npt.hpp" @@ -125,6 +127,11 @@ bool integrator_step_1(ParticleRange &particles) { // the Ermak-McCammon's Brownian Dynamics requires a single step // so, just skip here break; +#ifdef STOKESIAN_DYNAMICS + case INTEG_METHOD_SD: + stokesian_dynamics_step_1(particles); + break; +#endif // STOKESIAN_DYNAMICS default: throw std::runtime_error("Unknown value for integ_switch"); } @@ -150,6 +157,11 @@ void integrator_step_2(ParticleRange &particles) { // the Ermak-McCammon's Brownian Dynamics requires a single step brownian_dynamics_propagator(brownian, particles); break; +#ifdef STOKESIAN_DYNAMICS + case INTEG_METHOD_SD: + stokesian_dynamics_step_2(particles); + break; +#endif // STOKESIAN_DYNAMICS default: throw std::runtime_error("Unknown value for INTEG_SWITCH"); } @@ -202,7 +214,6 @@ int integrate(int n_steps, int reuse_forces) { #ifdef VALGRIND_INSTRUMENTATION CALLGRIND_START_INSTRUMENTATION; #endif - /* Integration loop */ ESPRESSO_PROFILER_CXX_MARK_LOOP_BEGIN(integration_loop, "Integration loop"); int integrated_steps = 0; @@ -280,6 +291,7 @@ int integrate(int n_steps, int reuse_forces) { } // for-loop over integration steps ESPRESSO_PROFILER_CXX_MARK_LOOP_END(integration_loop); + #ifdef VALGRIND_INSTRUMENTATION CALLGRIND_STOP_INSTRUMENTATION; #endif @@ -418,6 +430,16 @@ void integrate_set_bd() { mpi_bcast_parameter(FIELD_INTEG_SWITCH); } +int integrate_set_sd() { + if (box_geo.periodic(0) || box_geo.periodic(1) || box_geo.periodic(2)) { + runtimeErrorMsg() << "Stokesian Dynamics requires periodicity 0 0 0\n"; + return ES_ERROR; + } + integ_switch = INTEG_METHOD_SD; + mpi_bcast_parameter(FIELD_INTEG_SWITCH); + return ES_OK; +} + #ifdef NPT int integrate_set_npt_isotropic(double ext_pressure, double piston, bool xdir_rescale, bool ydir_rescale, diff --git a/src/core/integrate.hpp b/src/core/integrate.hpp index aee38748e44..6bd8e859819 100644 --- a/src/core/integrate.hpp +++ b/src/core/integrate.hpp @@ -31,6 +31,7 @@ #define INTEG_METHOD_NVT 1 #define INTEG_METHOD_STEEPEST_DESCENT 2 #define INTEG_METHOD_BD 3 +#define INTEG_METHOD_SD 7 /** Switch determining which integrator to use. */ extern int integ_switch; @@ -126,6 +127,9 @@ void integrate_set_nvt(); /** @brief Set the Brownian Dynamics integrator. */ void integrate_set_bd(); +/** @brief Set the Stokesian Dynamics integrator. */ +int integrate_set_sd(); + /** @brief Set the velocity Verlet integrator modified for the NpT ensemble * with isotropic rescaling. * diff --git a/src/core/integrators/stokesian_dynamics_inline.hpp b/src/core/integrators/stokesian_dynamics_inline.hpp new file mode 100644 index 00000000000..f44268942b7 --- /dev/null +++ b/src/core/integrators/stokesian_dynamics_inline.hpp @@ -0,0 +1,65 @@ +/* + * Copyright (C) 2010-2019 The ESPResSo project + * + * This file is part of ESPResSo. + * + * ESPResSo is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * ESPResSo is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ +#ifndef STOKESIAN_DYNAMICS_INLINE_HPP +#define STOKESIAN_DYNAMICS_INLINE_HPP + +#include "ParticleRange.hpp" +#include "config.hpp" +#include "particle_data.hpp" +#include "rotation.hpp" +#include "stokesian_dynamics/sd_interface.hpp" + +#include + +#ifdef STOKESIAN_DYNAMICS + +inline void +stokesian_dynamics_propagate_vel_pos(const ParticleRange &particles) { + auto const skin2 = Utils::sqr(0.5 * skin); + + // Compute new (translational and rotational) velocities + propagate_vel_pos_sd(particles); + + for (auto &p : particles) { + + // Translate particle + p.r.p += p.m.v * time_step; + + // Perform rotation + double norm = p.m.omega.norm(); + if (norm != 0) { + Utils::Vector3d omega_unit = (1 / norm) * p.m.omega; + local_rotate_particle(p, omega_unit, norm * time_step); + } + + // Verlet criterion check + if ((p.r.p - p.l.p_old).norm2() > skin2) + cell_structure.set_resort_particles(Cells::RESORT_LOCAL); + } +} + +inline void stokesian_dynamics_step_1(const ParticleRange &particles) { + stokesian_dynamics_propagate_vel_pos(particles); + sim_time += time_step; +} + +inline void stokesian_dynamics_step_2(const ParticleRange &particles) {} + +#endif // STOKESIAN_DYNAMICS +#endif diff --git a/src/core/serialization/SD_particle_data.hpp b/src/core/serialization/SD_particle_data.hpp new file mode 100644 index 00000000000..84fd2162839 --- /dev/null +++ b/src/core/serialization/SD_particle_data.hpp @@ -0,0 +1,56 @@ +/* + * Copyright (C) 2010-2020 The ESPResSo project + * + * This file is part of ESPResSo. + * + * ESPResSo is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * ESPResSo is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ +#ifndef SERIALIZATION_SD_PARTICLE_DATA_HPP +#define SERIALIZATION_SD_PARTICLE_DATA_HPP + +#include "stokesian_dynamics/sd_interface.hpp" + +#include +#include + +#if defined(STOKESIAN_DYNAMICS) or defined(STOKESIAN_DYNAMICS_GPU) + +BOOST_IS_BITWISE_SERIALIZABLE(SD_particle_data) + +namespace boost { +namespace serialization { +template +void load(Archive &ar, SD_particle_data &p, + const unsigned int /* file_version */) { + ar >> make_array(reinterpret_cast(&p), sizeof(SD_particle_data)); +} + +template +void save(Archive &ar, SD_particle_data const &p, + const unsigned int /* file_version */) { + /* Cruel but effective */ + ar << make_array(reinterpret_cast(&p), + sizeof(SD_particle_data)); +} + +template +void serialize(Archive &ar, SD_particle_data &p, + const unsigned int file_version) { + split_free(ar, p, file_version); +} +} // namespace serialization +} // namespace boost + +#endif +#endif diff --git a/src/core/stokesian_dynamics/CMakeLists.txt b/src/core/stokesian_dynamics/CMakeLists.txt new file mode 100644 index 00000000000..9292f026a36 --- /dev/null +++ b/src/core/stokesian_dynamics/CMakeLists.txt @@ -0,0 +1,25 @@ +# +# Copyright (C) 2019-2020 The ESPResSo project +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# + +target_include_directories( + EspressoCore PRIVATE ${stokesian_dynamics_SOURCE_DIR}/include + ${CMAKE_SOURCE_DIR}/src/core/stokesian_dynamics) + +target_sources(EspressoCore + PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/sd_interface.cpp) diff --git a/src/core/stokesian_dynamics/sd_interface.cpp b/src/core/stokesian_dynamics/sd_interface.cpp new file mode 100644 index 00000000000..bfb8ef08e31 --- /dev/null +++ b/src/core/stokesian_dynamics/sd_interface.cpp @@ -0,0 +1,308 @@ +/* + * Copyright (C) 2010-2020 The ESPResSo project + * + * This file is part of ESPResSo. + * + * ESPResSo is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * ESPResSo is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include "config.hpp" + +#if defined(STOKESIAN_DYNAMICS) || defined(STOKESIAN_DYNAMICS_GPU) + +#include +#include +#include + +#include + +#include "ParticleRange.hpp" +#include "communication.hpp" +#include "errorhandling.hpp" +#include "integrate.hpp" +#include "particle_data.hpp" +#include "rotation.hpp" +#include "thermostat.hpp" + +#include "serialization/SD_particle_data.hpp" + +#include +#include +#include + +#ifdef STOKESIAN_DYNAMICS +#include "stokesian_dynamics/sd_cpu.hpp" +#endif + +#ifdef STOKESIAN_DYNAMICS_GPU +#include "stokesian_dynamics/sd_gpu.hpp" +#endif + +#include "sd_interface.hpp" + +namespace { +double sd_viscosity = -1.0; + +enum { CPU, GPU, INVALID } device = INVALID; + +std::unordered_map radius_dict; + +double sd_kT = 0.0; + +std::size_t sd_seed = 0UL; + +int sd_flags = 0; + +/** Buffer that holds local particle data, and all particles on the master node + * used for sending particle data to master node. */ +std::vector parts_buffer{}; + +/** Buffer that holds the (translational and angular) velocities of the local + * particles on each node, used for returning results. */ +std::vector v_sd{}; + +} // namespace + +/** Pack selected properties of all local particles into a buffer. */ +void sd_gather_local_particles(ParticleRange const &parts) { + std::size_t n_part = parts.size(); + + parts_buffer.resize(n_part); + std::size_t i = 0; + + for (auto const &p : parts) { + parts_buffer[i].id = p.p.is_virtual ? -1 : p.p.identity; + parts_buffer[i].type = p.p.type; + + parts_buffer[i].pos[0] = p.r.p[0]; + parts_buffer[i].pos[1] = p.r.p[1]; + parts_buffer[i].pos[2] = p.r.p[2]; + + // Velocity is not needed for SD, thus omit while gathering. + // Particle orientation is also not needed for SD + // since all particles are assumed spherical. + + parts_buffer[i].ext_force[0] = p.f.f[0]; + parts_buffer[i].ext_force[1] = p.f.f[1]; + parts_buffer[i].ext_force[2] = p.f.f[2]; + + parts_buffer[i].ext_torque[0] = p.f.torque[0]; + parts_buffer[i].ext_torque[1] = p.f.torque[1]; + parts_buffer[i].ext_torque[2] = p.f.torque[2]; + + // radius_dict is not initialized on slave nodes -> need to assign radius + // later on master node + parts_buffer[i].r = 0; + + i++; + } +} + +/** Update translational and rotational velocities of all particles. */ +void sd_update_locally(ParticleRange const &parts) { + std::size_t i = 0; + + // Even though on the master node, the v_sd vector is larger than + // the (local) parts vector, this should still work. Because the local + // particles correspond to the first 6*n entries in the master's v_sd + // (which holds the velocities of ALL particles). + + for (auto &p : parts) { + // skip virtual particles + if (p.p.is_virtual) { + continue; + } + + // Copy velocities + p.m.v[0] = v_sd[6 * i]; + p.m.v[1] = v_sd[6 * i + 1]; + p.m.v[2] = v_sd[6 * i + 2]; + + p.m.omega[0] = v_sd[6 * i + 3]; + p.m.omega[1] = v_sd[6 * i + 4]; + p.m.omega[2] = v_sd[6 * i + 5]; + + i++; + } +} + +void set_sd_viscosity(double eta) { sd_viscosity = eta; } + +double get_sd_viscosity() { return sd_viscosity; } + +void set_sd_device(std::string const &dev) { + device = INVALID; +#ifdef STOKESIAN_DYNAMICS + if (dev == "cpu") { + device = CPU; + } +#endif +#ifdef STOKESIAN_DYNAMICS_GPU + if (dev == "gpu") { + device = GPU; + } +#endif +} + +std::string get_sd_device() { + switch (device) { + case CPU: + return "cpu"; + case GPU: + return "gpu"; + default: + return "invalid"; + } +} + +void set_sd_radius_dict(std::unordered_map const &x) { + radius_dict = x; +} + +std::unordered_map get_sd_radius_dict() { return radius_dict; } + +void set_sd_kT(double kT) { sd_kT = kT; } + +double get_sd_kT() { return sd_kT; } + +void set_sd_seed(std::size_t seed) { sd_seed = seed; } + +std::size_t get_sd_seed() { return sd_seed; } + +void set_sd_flags(int flg) { sd_flags = flg; } + +int get_sd_flags() { return sd_flags; } + +void propagate_vel_pos_sd(const ParticleRange &particles) { + std::size_t n_part_local = particles.size(); + + if (this_node == 0) { + if (integ_switch & INTEG_METHOD_SD) { + if (BOOST_UNLIKELY(sd_viscosity < 0.0)) { + runtimeErrorMsg() << "sd_viscosity has an invalid value: " + + std::to_string(sd_viscosity); + return; + } + + if (BOOST_UNLIKELY(sd_kT < 0.0)) { + runtimeErrorMsg() << "sd_kT has an invalid value: " + + std::to_string(sd_viscosity); + return; + } + + sd_gather_local_particles(particles); + Utils::Mpi::gather_buffer(parts_buffer, comm_cart, 0); + + std::size_t n_part = parts_buffer.size(); + + static std::vector id{}; + static std::vector x_host{}; + static std::vector f_host{}; + static std::vector a_host{}; + + // n_part not expected to change, but anyway ... + id.resize(n_part); + x_host.resize(6 * n_part); + f_host.resize(6 * n_part); + a_host.resize(n_part); + + std::size_t i = 0; + for (auto const &p : parts_buffer) { + id[i] = parts_buffer[i].id; + + x_host[6 * i + 0] = p.pos[0]; + x_host[6 * i + 1] = p.pos[1]; + x_host[6 * i + 2] = p.pos[2]; + // Actual orientation is not needed, just need default. + x_host[6 * i + 3] = 1; + x_host[6 * i + 4] = 0; + x_host[6 * i + 5] = 0; + + f_host[6 * i + 0] = p.ext_force[0]; + f_host[6 * i + 1] = p.ext_force[1]; + f_host[6 * i + 2] = p.ext_force[2]; + + f_host[6 * i + 3] = p.ext_torque[0]; + f_host[6 * i + 4] = p.ext_torque[1]; + f_host[6 * i + 5] = p.ext_torque[2]; + + double radius = radius_dict[p.type]; + + if (BOOST_UNLIKELY(radius < 0.0)) { + runtimeErrorMsg() + << "particle of type " + + std::to_string(int(parts_buffer[i].type)) + + " has an invalid radius: " + std::to_string(radius); + return; + } + + a_host[i] = radius; + + ++i; + } + + std::size_t offset = std::round(sim_time / time_step); + switch (device) { + +#ifdef STOKESIAN_DYNAMICS + case CPU: + v_sd = sd_cpu(x_host, f_host, a_host, n_part, sd_viscosity, + std::sqrt(sd_kT / time_step), offset, sd_seed, sd_flags); + break; +#endif + +#ifdef STOKESIAN_DYNAMICS_GPU + case GPU: + v_sd = sd_gpu(x_host, f_host, a_host, n_part, sd_viscosity, + std::sqrt(sd_kT / time_step), offset, sd_seed, sd_flags); + break; +#endif + + default: + runtimeErrorMsg() + << "Invalid device for Stokesian dynamics. Available devices:" +#ifdef STOKESIAN_DYNAMICS + " cpu" +#endif +#ifdef STOKESIAN_DYNAMICS_GPU + " gpu" +#endif + ; + return; + } + + Utils::Mpi::scatter_buffer( + v_sd.data(), static_cast(n_part_local * 6), comm_cart, 0); + sd_update_locally(particles); + } + + } else { // if (this_node == 0) + + if (thermo_switch & THERMO_SD) { + sd_gather_local_particles(particles); + Utils::Mpi::gather_buffer(parts_buffer, comm_cart, 0); + + v_sd.resize(n_part_local * 6); + + // now wait while master node is busy ... + + Utils::Mpi::scatter_buffer( + v_sd.data(), static_cast(n_part_local * 6), comm_cart, 0); + sd_update_locally(particles); + } + + } // if (this_node == 0) {...} else +} + +#endif // STOKESIAN_DYNAMICS diff --git a/src/core/stokesian_dynamics/sd_interface.hpp b/src/core/stokesian_dynamics/sd_interface.hpp new file mode 100644 index 00000000000..e36a0efe520 --- /dev/null +++ b/src/core/stokesian_dynamics/sd_interface.hpp @@ -0,0 +1,89 @@ +/* + * Copyright (C) 2010-2020 The ESPResSo project + * + * This file is part of ESPResSo. + * + * ESPResSo is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * ESPResSo is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +/** @file + * See @cite durlofsky87a for the Stokesian dynamics method used here. + * See @cite banchio03a and @cite brady88a for the thermalization method. + */ + +#ifndef STOKESIAN_DYNAMICS_INTERFACE_H +#define STOKESIAN_DYNAMICS_INTERFACE_H + +#include "ParticleRange.hpp" +#include "config.hpp" + +#include +#include +#include + +#ifdef STOKESIAN_DYNAMICS + +/* type for particle data transfer between nodes */ +struct SD_particle_data { + int on_node = -1; + int id = -1; + int type = 0; + + /* particle radius */ + double r = -1; + + /* particle position */ + Utils::Vector3d pos = {0., 0., 0.}; + + /* particle velocity */ + Utils::Vector3d vel = {0., 0., 0.}; + + /* particle rotational velocity */ + Utils::Vector3d omega = {0., 0., 0.}; + + /* external force */ + Utils::Vector3d ext_force = {0.0, 0.0, 0.0}; + + /* external torque */ + Utils::Vector3d ext_torque = {0.0, 0.0, 0.0}; +}; + +void set_sd_viscosity(double eta); +double get_sd_viscosity(); + +void set_sd_device(std::string const &dev); +std::string get_sd_device(); + +void set_sd_radius_dict(std::unordered_map const &x); +std::unordered_map get_sd_radius_dict(); + +void set_sd_kT(double kT); +double get_sd_kT(); + +void set_sd_flags(int flg); +int get_sd_flags(); + +void set_sd_seed(std::size_t seed); +std::size_t get_sd_seed(); + +/** Takes the forces and torques on all particles and computes their + * velocities. Acts globally on particles on all nodes; i.e. particle data + * is gathered from all nodes and their velocities and angular velocities are + * set according to the Stokesian Dynamics method. + */ +void propagate_vel_pos_sd(const ParticleRange &particles); + +#endif // STOKESIAN_DYNAMICS + +#endif // STOKESIAN_DYNAMICS_INTERFACE_H diff --git a/src/core/thermostat.hpp b/src/core/thermostat.hpp index cc7115a4e83..f87b600ede1 100644 --- a/src/core/thermostat.hpp +++ b/src/core/thermostat.hpp @@ -46,6 +46,7 @@ #define THERMO_NPT_ISO 4 #define THERMO_LB 8 #define THERMO_BROWNIAN 16 +#define THERMO_SD 32 /*@}*/ namespace Thermostat { diff --git a/src/core/unit_tests/reaction_ensemble_utils_test.cpp b/src/core/unit_tests/reaction_ensemble_utils_test.cpp index 29b04c0db03..90664cc98bd 100644 --- a/src/core/unit_tests/reaction_ensemble_utils_test.cpp +++ b/src/core/unit_tests/reaction_ensemble_utils_test.cpp @@ -39,8 +39,8 @@ BOOST_AUTO_TEST_CASE(average_list_of_allowed_entries_test) { using namespace ReactionEnsemble; constexpr double tol = 100 * std::numeric_limits::epsilon(); - BOOST_CHECK_CLOSE( - average_list_of_allowed_entries(std::vector{1, 2}), 1.5, tol); + BOOST_CHECK_CLOSE(average_list_of_allowed_entries(std::vector{1, 2}), + 1.5, tol); BOOST_CHECK_CLOSE( average_list_of_allowed_entries(std::vector{1, 2, -2}), 1.5, tol); BOOST_CHECK_CLOSE( diff --git a/src/python/espressomd/globals.pyx b/src/python/espressomd/globals.pyx index 719cf68761c..20ca445af16 100644 --- a/src/python/espressomd/globals.pyx +++ b/src/python/espressomd/globals.pyx @@ -25,7 +25,7 @@ from .globals cimport timing_samples from .globals cimport forcecap_set from .globals cimport forcecap_get from .utils import array_locked -from .utils cimport Vector3d, make_array_locked +from .utils cimport Vector3d, make_array_locked, handle_errors cdef class Globals: property box_l: @@ -68,6 +68,7 @@ cdef class Globals: box_geo.set_periodic(i, _periodic[i]) mpi_bcast_parameter(FIELD_PERIODIC) + handle_errors("Error while assigning system periodicity") def __get__(self): periodicity = np.zeros(3) diff --git a/src/python/espressomd/integrate.pxd b/src/python/espressomd/integrate.pxd index 60a6577adfb..6fe47ebeaa2 100644 --- a/src/python/espressomd/integrate.pxd +++ b/src/python/espressomd/integrate.pxd @@ -17,6 +17,8 @@ # along with this program. If not, see . # from libcpp cimport bool as cbool +from libcpp.string cimport string +from libcpp.unordered_map cimport unordered_map include "myconfig.pxi" @@ -25,6 +27,7 @@ cdef extern from "config.hpp": cdef extern from "integrate.hpp" nogil: cdef int python_integrate(int n_steps, cbool recalc_forces, int reuse_forces) + cdef void integrate_set_sd() cdef void integrate_set_nvt() cdef int integrate_set_steepest_descent(const double f_max, const double gamma, const int max_steps, const double max_displacement) @@ -38,6 +41,28 @@ IF NPT: cbool xdir_rescale, cbool ydir_rescale, cbool zdir_rescale, cbool cubic_box) +cdef extern from "stokesian_dynamics/sd_interface.hpp": + IF(STOKESIAN_DYNAMICS or STOKESIAN_DYNAMICS_GPU): + void set_sd_viscosity(double eta) + double get_sd_viscosity() + + void set_sd_device(const string & dev) + string get_sd_device() + + void set_sd_radius_dict(const unordered_map[int, double] & radius_dict) + unordered_map[int, double] get_sd_radius_dict() + + void set_sd_flags(int flg) + int get_sd_flags() + +IF(STOKESIAN_DYNAMICS or STOKESIAN_DYNAMICS_GPU): + cpdef enum flags: + NONE = 0, + SELF_MOBILITY = 1 << 0, + PAIR_MOBILITY = 1 << 1, + LUBRICATION = 1 << 2, + FTS = 1 << 3 + cdef inline int _integrate(int nSteps, cbool recalc_forces, int reuse_forces): with nogil: return python_integrate(nSteps, recalc_forces, reuse_forces) diff --git a/src/python/espressomd/integrate.pyx b/src/python/espressomd/integrate.pyx index 044a4640dd1..9b698850f45 100644 --- a/src/python/espressomd/integrate.pyx +++ b/src/python/espressomd/integrate.pyx @@ -19,6 +19,7 @@ from cpython.exc cimport PyErr_CheckSignals, PyErr_SetInterrupt include "myconfig.pxi" from .utils cimport handle_errors, check_type_or_throw_except +from .utils import to_char_pointer from . cimport integrate cdef class IntegratorHandle: @@ -100,6 +101,13 @@ cdef class IntegratorHandle: """ self._integrator = BrownianDynamics() + def set_stokesian_dynamics(self, *args, **kwargs): + """ + Set the integration method to Stokesian Dynamics (:class:`StokesianDynamics`). + + """ + self._integrator = StokesianDynamics(*args, **kwargs) + cdef class Integrator: """ @@ -405,3 +413,105 @@ cdef class BrownianDynamics(Integrator): def _set_params_in_es_core(self): integrate_set_bd() + + +IF(STOKESIAN_DYNAMICS or STOKESIAN_DYNAMICS_GPU): + cdef class StokesianDynamics(Integrator): + """ + Stokesian Dynamics integrator. + + Parameters + ---------- + viscosity : :obj:`float` + Bulk viscosity. + radii : :obj:`dict` + Dictionary that maps particle types to radii. + device : :obj:`str`, optional, \{'cpu', 'gpu'\} + Device to execute on. + approximation_method : :obj:`str`, optional, \{'ft', 'fts'\} + Chooses the method of the mobility approximation. + ``'fts'`` is more accurate. Default is ``'fts'``. + self_mobility : :obj:`bool`, optional + Switches off or on the mobility terms for single particles. Default + is ``True``. + pair_mobility : :obj:`bool`, optional + Switches off or on the hydrodynamic interactions between particles. + Default is ``True``. + + """ + + def default_params(self): + IF STOKESIAN_DYNAMICS: + sd_device_str = "cpu" + ELIF STOKESIAN_DYNAMICS_GPU: + sd_device_str = "gpu" + return {"lubrication": False, "approximation_method": "fts", + "self_mobility": True, "pair_mobility": True, + "device": sd_device_str} + + def valid_keys(self): + """All parameters that can be set. + + """ + return {"radii", "viscosity", "device", "lubrication", + "approximation_method", "self_mobility", "pair_mobility"} + + def required_keys(self): + """Parameters that have to be set. + + """ + return {"radii", "viscosity"} + + def validate_params(self): + check_type_or_throw_except( + self._params["viscosity"], 1, float, + "viscosity must be a number") + check_type_or_throw_except( + self._params["device"], 1, str, + "device must be a string") + check_type_or_throw_except( + self._params["radii"], 1, dict, + "radii must be a dictionary") + check_type_or_throw_except( + self._params["lubrication"], 1, bool, + "lubrication must be a bool") + if self._params["lubrication"]: + raise NotImplementedError( + "Stokesian Dynamics lubrication is not available yet") + check_type_or_throw_except( + self._params["approximation_method"], 1, str, + "approximation_method must be a string") + if self._params["approximation_method"].lower() not in { + "ft", "fts"}: + raise ValueError( + "approximation_method must be either 'ft' or 'fts'") + check_type_or_throw_except( + self._params["self_mobility"], 1, bool, + "self_mobility must be a bool") + check_type_or_throw_except( + self._params["pair_mobility"], 1, bool, + "pair_mobility must be a bool") + + def _set_params_in_es_core(self): + integrate_set_sd() + set_sd_radius_dict(self._params["radii"]) + set_sd_device(to_char_pointer(self._params["device"].lower())) + set_sd_viscosity(self._params["viscosity"]) + fl = flags.NONE + if self._params["lubrication"]: + fl = fl | flags.LUBRICATION + if self._params["approximation_method"].lower() == "fts": + fl = fl | flags.FTS + if self._params["self_mobility"]: + fl = fl | flags.SELF_MOBILITY + if self._params["pair_mobility"]: + fl = fl | flags.PAIR_MOBILITY + set_sd_flags(fl) + + handle_errors( + "Encountered error while setting integration method to SD") + +ELSE: + cdef class StokesianDynamics(Integrator): + def __init__(self, *args, **kwargs): + raise Exception("Stokesian Dynamics not compiled in.") diff --git a/src/python/espressomd/thermostat.pxd b/src/python/espressomd/thermostat.pxd index 707cd2de124..622cda0f333 100644 --- a/src/python/espressomd/thermostat.pxd +++ b/src/python/espressomd/thermostat.pxd @@ -32,6 +32,7 @@ cdef extern from "thermostat.hpp": int THERMO_NPT_ISO int THERMO_DPD int THERMO_BROWNIAN + int THERMO_SD IF PARTICLE_ANISOTROPY: ctypedef struct langevin_thermostat_struct "LangevinThermostat": @@ -69,6 +70,14 @@ cdef extern from "thermostat.hpp": IF DPD: stdint.uint64_t dpd_get_rng_state() +cdef extern from "stokesian_dynamics/sd_interface.hpp": + IF(STOKESIAN_DYNAMICS or STOKESIAN_DYNAMICS_GPU): + void set_sd_kT(double kT) + double get_sd_kT() + + void set_sd_seed(size_t seed) + size_t get_sd_seed() + cdef extern from "script_interface/Globals.hpp": # links intern C-struct with python object cdef extern langevin_thermostat_struct langevin diff --git a/src/python/espressomd/thermostat.pyx b/src/python/espressomd/thermostat.pyx index a89bf3994ac..2636ac0255e 100644 --- a/src/python/espressomd/thermostat.pyx +++ b/src/python/espressomd/thermostat.pyx @@ -124,6 +124,8 @@ cdef class Thermostat: gamma_rotation=thmst["gamma_rotation"], act_on_virtual=thmst["act_on_virtual"], seed=thmst["seed"]) + if thmst["type"] == "SD": + self.set_stokesian(kT=thmst["kT"], seed=thmst["seed"]) def get_ts(self): return thermo_switch @@ -205,6 +207,13 @@ cdef class Thermostat: dpd_dict["kT"] = temperature dpd_dict["seed"] = int(dpd_get_rng_state()) thermo_list.append(dpd_dict) + if (thermo_switch & THERMO_SD): + IF STOKESIAN_DYNAMICS: + sd_dict = {} + sd_dict["type"] = "SD" + sd_dict["kT"] = get_sd_kT() + sd_dict["seed"] = get_sd_seed() + thermo_list.append(sd_dict) return thermo_list def turn_off(self): @@ -698,3 +707,35 @@ cdef class Thermostat: mpi_bcast_parameter(FIELD_THERMO_SWITCH) mpi_bcast_parameter(FIELD_TEMPERATURE) + + IF STOKESIAN_DYNAMICS: + def set_stokesian(self, kT=None, seed=None): + """ + Sets the SD thermostat with required parameters. + + This thermostat requires the feature ``STOKESIAN_DYNAMICS`` or + ``STOKESIAN_DYNAMICS_GPU``. + + Parameters + ---------- + kT : :obj:`float`, optional + Temperature + seed : :obj:`int`, optional + Seed for the random number generator + + """ + + if (kT is None) or (kT == 0): + set_sd_kT(0.0) + else: + utils.check_type_or_throw_except( + kT, 1, float, "kT must be a float") + set_sd_kT(kT) + + utils.check_type_or_throw_except( + seed, 1, int, "seed must be an integer") + set_sd_seed(seed) + + global thermo_switch + thermo_switch = (thermo_switch | THERMO_SD) + mpi_bcast_parameter(FIELD_THERMO_SWITCH) diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index ee1d8d0646a..f98094a6ede 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -56,7 +56,7 @@ endfunction(PYTHON_TEST) # Separate features with hyphens, use a period to add an optional flag. foreach( TEST_COMBINATION - lb.cpu-p3m.cpu-lj-therm.lb;lb.gpu-p3m.cpu-lj-therm.lb;ek.gpu;lb.off-therm.npt-int.npt-minimization;lb.off-therm.langevin-int.nvt;lb.off-therm.dpd-int.sd;lb.off-therm.bd-int.bd + lb.cpu-p3m.cpu-lj-therm.lb;lb.gpu-p3m.cpu-lj-therm.lb;ek.gpu;lb.off-therm.npt-int.npt-minimization;lb.off-therm.langevin-int.nvt;lb.off-therm.dpd-int.sd;lb.off-therm.bd-int.bd;lb.off-therm.sdm-int.sdm.cpu;lb.off-therm.sdm-int.sdm.gpu ) if(${TEST_COMBINATION} MATCHES "\\.gpu") set(TEST_LABELS "gpu") @@ -200,9 +200,11 @@ python_test(FILE gpu_availability.py MAX_NUM_PROC 1 LABELS gpu) python_test(FILE features.py MAX_NUM_PROC 1) python_test(FILE galilei.py MAX_NUM_PROC 32) python_test(FILE time_series.py MAX_NUM_PROC 1) -python_test(FILE linear_momentum.py) +python_test(FILE linear_momentum.py MAX_NUM_PROC 4) python_test(FILE linear_momentum_lb.py MAX_NUM_PROC 2 LABELS gpu) python_test(FILE mmm1d.py MAX_NUM_PROC 2 LABELS gpu) +python_test(FILE stokesian_dynamics_cpu.py MAX_NUM_PROC 2) +python_test(FILE stokesian_dynamics_gpu.py MAX_NUM_PROC 4 LABELS gpu) python_test(FILE elc.py MAX_NUM_PROC 2) python_test(FILE elc_vs_analytic.py MAX_NUM_PROC 2) python_test(FILE rotation.py MAX_NUM_PROC 1) @@ -225,6 +227,9 @@ add_custom_target( ${CMAKE_CURRENT_BINARY_DIR} COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/ek_common.py ${CMAKE_CURRENT_BINARY_DIR} + COMMAND + ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/stokesian_dynamics.py + ${CMAKE_CURRENT_BINARY_DIR} COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/ek_eof_one_species_base.py diff --git a/testsuite/python/data/dancing.txt b/testsuite/python/data/dancing.txt new file mode 100644 index 00000000000..e4628ebf48f --- /dev/null +++ b/testsuite/python/data/dancing.txt @@ -0,0 +1,125 @@ +# F-T-S trajectory digitized from Figure 5b in Durlofsky, Brady, Bossis +# (1987). Dynamic simulation of hydrodynamically interacting particles. +# Journal of Fluid Mechanics, 180, 21-49. doi:10.1017/S002211208700171X +-5.0 0.0 0.0 0.0 7.0 0.0 +-4.99668 -0.713 -0.01829 -0.487 6.94144 -0.271 +-4.98734 -6.788 -0.02043 -6.451 6.96641 -5.903 +-4.96921 -12.86 -0.02558 -10.758 6.8727 -11.347 +-4.90709 -19.304 -0.03301 -20.568 6.77196 -16.189 +-4.82728 -26.884 -0.02689 -27.026 6.67141 -23.141 +-4.6861 -32.553 -0.0403 -31.814 6.51147 -29.351 +-4.54493 -38.222 -0.04236 -39.709 6.24687 -35.582 +-4.36852 -44.644 -0.00852 -46.52 5.99622 -41.659 +-4.17456 -50.682 0.01422 -52.975 5.7281 -47.588 +-3.95402 -58.994 0.07287 -58.586 5.41122 -54.431 +-3.71617 -64.265 0.104 -66.752 5.07323 -59.32 +-3.52218 -70.683 0.23208 -73.373 4.71101 -66.322 +-3.2843 -76.334 0.30747 -80.764 4.33477 -72.574 +-3.08145 -83.51 0.43989 -86.815 3.98316 -81.233 +-2.95778 -89.942 0.57249 -94.954 3.63823 -86.424 +-2.8606 -95.239 0.74897 -102.135 3.26553 -93.278 +-2.78066 -104.338 0.91214 -107.99 2.94153 -98.616 +-2.77128 -110.792 1.11486 -113.647 2.61085 -107.119 +-2.73554 -117.241 1.33543 -122.339 2.29731 -112.304 +-2.72613 -124.075 1.58657 -128.746 1.98406 -120.654 +-2.72545 -131.67 1.82438 -133.637 1.73691 -127.032 +-2.76869 -139.654 2.11099 -143.455 1.4199 -132.368 +-2.80344 -144.218 2.42352 -148.141 1.17984 -139.95 +-2.85537 -153.342 2.73175 -153.968 0.92575 -146.63 +-2.97797 -158.682 3.06642 -160.549 0.71708 -154.055 +-3.09168 -165.16 3.40558 -168.078 0.57461 -160.413 +-3.24062 -170.886 3.8016 -172.559 0.40082 -167.68 +-3.4247 -176.997 4.18904 -179.319 0.25487 -174.189 +-3.55606 -182.719 4.58096 -187.028 0.14042 -181.897 +-3.79271 -190.74 4.98141 -191.888 0.06416 -187.338 +-3.96795 -197.23 5.31609 -198.469 0.05085 -194.574 +-4.11683 -203.715 5.66836 -205.236 0.04098 -201.207 +-4.20433 -208.289 5.9678 -211.064 0.12173 -206.768 +-4.27395 -216.277 6.24972 -217.655 0.18532 -215.497 +-4.2736 -220.075 6.54027 -222.346 0.28694 -220.45 +-4.27296 -227.29 6.76055 -227.81 0.40972 -227.811 +-4.25479 -233.743 6.90168 -232.909 0.48357 -234.277 +-4.17502 -240.943 7.04732 -239.337 0.63075 -241.483 +-4.10424 -245.866 7.19739 -246.144 0.75692 -247.788 +-3.95417 -252.673 7.24187 -252.021 0.88652 -253.49 +-3.74263 -258.708 7.26443 -258.473 0.99859 -258.14 +-3.56612 -266.269 7.19456 -263.613 1.00915 -259.193 +-3.38095 -272.309 7.11152 -268.946 1.13539 -266.252 +-3.14307 -277.959 6.98464 -275.616 1.19164 -271.064 +-2.83482 -283.976 6.76989 -282.684 1.30061 -280.085 +-2.55295 -289.997 6.54177 -287.855 1.33956 -286.407 +-2.22712 -296.01 6.20834 -295.135 1.44128 -292.416 +-1.93656 -300.891 5.84403 -301.092 1.45917 -297.386 +-1.65444 -309.571 5.52372 -307.61 1.53676 -306.564 +-1.43408 -315.984 5.20323 -312.039 1.58247 -310.473 +-1.24005 -322.782 4.86109 -320.27 1.64259 -319.504 +-1.09891 -328.071 4.56718 -326.973 1.69221 -328.085 +-1.00155 -335.268 4.21171 -333.498 1.7311 -333.804 +-0.94812 -342.853 3.84755 -341.164 1.84686 -340.715 +-0.92989 -350.065 3.5536 -347.487 1.9244 -349.441 +-0.97323 -356.909 3.22455 -354.576 2.08887 -354.834 +-0.99016 -364.127 2.93498 -360.708 2.20809 -361.442 +-1.09512 -370.224 2.64551 -367.98 2.43208 -369.236 +-1.12963 -377.446 2.37366 -375.628 2.63853 -375.978 +-1.20821 -383.538 2.12367 -381.942 2.90775 -382.255 +-1.2691 -390.765 1.91341 -390.148 3.20145 -389.281 +-1.3123 -399.128 1.80408 -396.435 3.57531 -395.237 +-1.3381 -405.589 1.75649 -404.989 3.92491 -403.157 +-1.27601 -411.653 1.84504 -411.997 4.41045 -409.543 +-1.10832 -418.836 2.05227 -418.983 4.89598 -415.779 +-0.84379 -427.519 2.45288 -425.551 5.39202 -422.615 +-0.42119 -434.273 2.92399 -434.004 5.94386 -429.14 +0.04537 -441.018 3.48293 -441.491 6.49245 -438.227 +0.55576 -446.236 4.20446 -448.377 7.03008 -441.89 +1.08382 -452.59 4.82041 -454.334 7.54701 -448.27 +1.61196 -459.703 5.40148 -463.525 8.05703 -455.255 +2.14878 -465.675 5.82411 -470.469 8.55671 -463.748 +2.65038 -470.894 6.19403 -477.992 8.93058 -469.855 +3.10828 -479.16 6.44075 -484.211 9.23483 -477.783 +3.55709 -484.01 6.5294 -492.358 9.42714 -482.418 +3.94463 -491.91 6.4949 -499.77 9.47669 -490.245 +4.34086 -498.669 6.39002 -506.626 9.46054 -504.566 +4.63156 -505.068 6.22348 -512.165 9.37751 -512.419 +4.93971 -509.946 6.00879 -519.992 9.31517 -517.555 +5.13387 -518.263 5.81163 -527.245 9.15885 -525.12 +5.31914 -525.442 5.5661 -534.319 9.06867 -531.016 +5.48672 -531.486 5.32066 -542.341 8.95083 -539.78 +5.60312 -537.658 5.02225 -547.906 8.85011 -544.923 +5.6895 -543.622 4.77684 -556.308 8.68331 -552.34 +5.7635 -550.485 4.52669 -560.914 8.55833 -559.297 +5.83875 -555.794 4.28115 -567.797 8.46127 -566.248 +5.85595 -562.25 4.03131 -575.82 8.3676 -572.144 +5.90105 -571.093 3.73739 -582.523 8.27406 -579.547 +5.89598 -576.357 3.52696 -588.83 8.16984 -584.541 +5.89652 -582.338 3.33401 -593.994 8.09043 -593.749 +5.89717 -589.514 3.1367 -599.539 8.03155 -598.583 +5.88109 -595.498 2.92653 -608.694 7.95536 -604.928 +5.89271 -601.477 2.73357 -613.858 7.89312 -611.119 +5.86555 -607.463 2.58029 -620.154 7.84145 -618.514 +5.8384 -613.449 2.42707 -627.209 7.79681 -626.51 +5.83336 -618.953 2.26923 -631.797 7.79035 -632.238 +5.80075 -625.897 2.12917 -638.47 7.75944 -637.669 +5.77909 -631.403 2.02867 -645.135 7.77052 -644.449 +5.7353 -637.153 1.84456 -650.867 7.77796 -649.572 +5.69149 -642.664 1.77029 -656.008 7.77162 -656.656 +5.66458 -651.282 1.61702 -662.494 7.80708 -662.979 +5.67063 -657.022 1.49023 -670.114 7.80755 -668.254 +5.61581 -663.253 1.35904 -677.734 7.83945 -673.824 +5.59983 -670.314 1.30229 -682.113 7.87851 -681.352 +5.54213 -675.349 1.21042 -687.067 7.93472 -685.711 +5.50668 -681.456 1.08801 -694.496 7.99107 -691.578 +5.47959 -688.16 0.97429 -700.785 8.051 -698.348 +5.44966 -694.027 0.90881 -705.924 8.06904 -704.976 +5.40864 -699.777 0.7995 -712.401 8.15679 -710.837 +5.37055 -707.32 0.70767 -717.736 8.22705 -716.248 +5.31297 -713.671 0.64245 -725.723 8.31842 -723.465 +5.30231 -718.338 0.50232 -731.637 8.38871 -729.178 +5.31399 -725.034 0.41486 -736.78 8.48694 -735.187 +5.27302 -731.382 0.35396 -743.818 8.58866 -741.196 +5.23482 -737.609 0.20513 -750.872 8.68342 -747.357 +5.21315 -742.996 0.15274 -754.87 8.80963 -754.114 +5.14998 -748.87 0.02588 -761.73 8.9148 -759.67 +5.15612 -755.448 -0.03085 -766.298 9.02002 -765.828 +5.11519 -762.274 -0.12256 -772.962 9.11844 -773.948 +5.11009 -767.179 -0.30214 -780.212 9.24798 -778.896 +5.09677 -772.923 -0.37187 -786.871 9.36707 -784.148 diff --git a/testsuite/python/save_checkpoint.py b/testsuite/python/save_checkpoint.py index 78ffe793041..1260a9d0715 100644 --- a/testsuite/python/save_checkpoint.py +++ b/testsuite/python/save_checkpoint.py @@ -153,6 +153,9 @@ system.thermostat.set_npt(kT=1.0, gamma0=2.0, gammav=0.1, seed=42) elif 'THERM.DPD' in modes and has_features('DPD'): system.thermostat.set_dpd(kT=1.0, seed=42) + elif 'THERM.SDM' in modes and (has_features('STOKESIAN_DYNAMICS') or has_features('STOKESIAN_DYNAMICS_GPU')): + system.periodicity = [0, 0, 0] + system.thermostat.set_stokesian(kT=1.0, seed=42) # set integrator if 'INT.NPT' in modes and has_features('NPT'): system.integrator.set_isotropic_npt(ext_pressure=2.0, piston=0.01, @@ -164,6 +167,16 @@ system.integrator.set_nvt() elif 'INT.BD' in modes: system.integrator.set_brownian_dynamics() + elif 'INT.SDM.CPU' in modes and has_features('STOKESIAN_DYNAMICS'): + system.periodicity = [0, 0, 0] + system.integrator.set_stokesian_dynamics( + approximation_method='ft', device='cpu', viscosity=0.5, + radii={0: 1.5}, pair_mobility=False, self_mobility=True) + elif 'INT.SDM.GPU' in modes and has_features('STOKESIAN_DYNAMICS_GPU') and espressomd.gpu_available(): + system.periodicity = [0, 0, 0] + system.integrator.set_stokesian_dynamics( + approximation_method='fts', device='gpu', viscosity=2.0, + radii={0: 1.0}, pair_mobility=True, self_mobility=False) # set minimization if 'MINIMIZATION' in modes: steepest_descent(system, f_max=1, gamma=10, max_steps=0, diff --git a/testsuite/python/stokesian_dynamics.py b/testsuite/python/stokesian_dynamics.py new file mode 100644 index 00000000000..3f710126421 --- /dev/null +++ b/testsuite/python/stokesian_dynamics.py @@ -0,0 +1,203 @@ +# +# Copyright (C) 2019-2020 The ESPResSo project +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +import espressomd +import espressomd.observables +import espressomd.accumulators +import espressomd.constraints +import numpy as np +import unittest as ut +from tests_common import abspath + +s = espressomd.System(box_l=[1.0, 1.0, 1.0]) + + +class StokesianDynamicsSetupTest(ut.TestCase): + system = s + device = 'none' + + def setUp(self): + self.system.box_l = [10] * 3 + + self.system.time_step = 1.0 + self.system.cell_system.skin = 0.4 + + # unset SD integrator so we can test whether set_stokesian_dynamics() + # fails. set_nvt() is the only way to ensure that integ_switch is + # set to a different value than INTEG_METHOD_SD + self.system.integrator.set_nvt() + + def pbc_checks(self): + self.system.periodicity = [0, 0, 1] + with self.assertRaises(Exception): + self.system.integrator.set_stokesian_dynamics( + viscosity=1.0, device=self.device, radii={0: 1.0}) + + self.system.periodicity = [0, 0, 0] + self.system.integrator.set_stokesian_dynamics( + viscosity=1.0, device=self.device, radii={0: 1.0}) + + with self.assertRaises(Exception): + self.system.periodicity = [0, 1, 0] + + +class StokesianDynamicsTest(ut.TestCase): + system = s + device = 'none' + + # Digitized reference data of Figure 5b from + # Durlofsky et al., J. Fluid Mech. 180, 21 (1987) + # https://doi.org/10.1017/S002211208700171X + data = np.loadtxt(abspath('data/dancing.txt')) + + def setUp(self): + self.system.box_l = [10] * 3 + self.system.periodicity = [0, 0, 0] + self.system.cell_system.skin = 0.4 + + def falling_spheres(self, time_step, l_factor, t_factor, + sd_method='fts', sd_short=False): + self.system.time_step = time_step + self.system.part.add(pos=[-5 * l_factor, 0, 0], rotation=[1, 1, 1]) + self.system.part.add(pos=[0 * l_factor, 0, 0], rotation=[1, 1, 1]) + self.system.part.add(pos=[7 * l_factor, 0, 0], rotation=[1, 1, 1]) + + self.system.integrator.set_stokesian_dynamics( + viscosity=1.0 / (t_factor * l_factor), + device=self.device, radii={0: 1.0 * l_factor}, + approximation_method=sd_method) + + gravity = espressomd.constraints.Gravity( + g=[0, -1.0 * l_factor / (t_factor**2), 0]) + self.system.constraints.add(gravity) + self.system.time_step = 1.0 * t_factor + + obs = espressomd.observables.ParticlePositions(ids=(0, 1, 2)) + acc = espressomd.accumulators.TimeSeries(obs=obs, delta_N=1) + self.system.auto_update_accumulators.add(acc) + acc.update() + + if sd_short: + y_min = -80 + intsteps = 1300 + elif sd_method == 'fts': + y_min = -555 + intsteps = 8000 + else: + y_min = -200 + intsteps = 3000 + intsteps = int(intsteps * t_factor / self.system.time_step) + + self.system.integrator.run(intsteps) + + simul = acc.time_series()[:, :, 0:2] + paper = self.data.reshape([-1, 3, 2]) + + for pid in range(3): + dist = [] + # the simulated trajectory is oversampled by a ratio of + # (90/t_factor):1 compared to the published trajectory + for desired in paper[:, pid] * l_factor: + if desired[1] < y_min * l_factor: + break + # find the closest point in the simulated trajectory + idx = np.abs(simul[:, pid, 1] - desired[1]).argmin() + actual = simul[idx, pid] + dist.append(np.linalg.norm(actual - desired)) + self.assertLess(idx, intsteps, msg='Insufficient sampling') + np.testing.assert_allclose(dist, 0, rtol=0, atol=0.5 * l_factor) + + def tearDown(self): + self.system.constraints.clear() + self.system.part.clear() + + +class StokesianDiffusionTest(ut.TestCase): + system = s + device = 'none' + + kT = 1e-4 + R = 1.5 + eta = 2.4 + + def setUp(self): + self.system.box_l = [10] * 3 + self.system.periodicity = [0, 0, 0] + + self.system.time_step = 1.0 + self.system.cell_system.skin = 0.4 + + self.system.part.add(pos=[0, 0, 0], rotation=[1, 1, 1]) + + def check(self): + self.system.integrator.set_stokesian_dynamics( + viscosity=self.eta, device=self.device, radii={0: self.R}) + self.system.thermostat.set_stokesian(kT=self.kT, seed=42) + + intsteps = int(100000 / self.system.time_step) + pos = np.empty([intsteps + 1, 3]) + orientation = np.empty((intsteps + 1, 3)) + + pos[0, :] = self.system.part[0].pos + orientation[0, :] = self.system.part[0].director + for i in range(intsteps): + self.system.integrator.run(1) + pos[i + 1, :] = self.system.part[0].pos + orientation[i + 1, :] = self.system.part[0].director + + # translational diffusion coefficient + D_expected = self.kT / (6 * np.pi * self.eta * self.R) + + # NOTE on steps_per_slice: + # The shorter these trajectories are, the more trajectories we get + # and the more accurate the result is. + # However since we want to test that diffusion works for + # "long" trajectories we won't go too low. + n_steps_per_slice = 200 + n_slices = int(intsteps / n_steps_per_slice) + + squared_displacement_per_slice = np.empty(n_slices) + for i in range(n_slices): + squared_displacement_per_slice[i] = np.linalg.norm( + pos[i * n_steps_per_slice] - pos[(i + 1) * n_steps_per_slice])**2 + D_measured = np.mean(squared_displacement_per_slice) / \ + (6 * n_steps_per_slice * self.system.time_step) + self.assertAlmostEqual(D_expected, D_measured, delta=D_expected * 0.1) + + # rotational diffusion coefficient + Dr_expected = self.kT / (8 * np.pi * self.eta * self.R**3) + n_steps_per_slice = 200 + n_slices = int(intsteps / n_steps_per_slice) + + # This is the displacement equivalent to measure rotational diffusion + costheta_per_slice = np.empty(n_slices) + + for i in range(n_slices): + costheta_per_slice[i] = np.dot(orientation[i * n_steps_per_slice, :], + orientation[(i + 1) * n_steps_per_slice, :]) + Dr_measured = -np.log(np.mean(costheta_per_slice)) / \ + (2 * n_steps_per_slice * self.system.time_step) + self.assertAlmostEqual( + Dr_expected, + Dr_measured, + delta=Dr_expected * 0.1) + + def tearDown(self): + self.system.constraints.clear() + self.system.part.clear() + self.system.thermostat.set_stokesian(kT=0) diff --git a/testsuite/python/stokesian_dynamics_cpu.py b/testsuite/python/stokesian_dynamics_cpu.py new file mode 100644 index 00000000000..4959da98631 --- /dev/null +++ b/testsuite/python/stokesian_dynamics_cpu.py @@ -0,0 +1,58 @@ +# +# Copyright (C) 2019-2020 The ESPResSo project +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +import unittest as ut +import unittest_decorators as utx +import stokesian_dynamics as sd + + +@utx.skipIfMissingFeatures(["STOKESIAN_DYNAMICS"]) +class StokesianDynamicsSetupTest(sd.StokesianDynamicsSetupTest): + device = 'cpu' + + def test_pbc_checks(self): + self.pbc_checks() + + +@utx.skipIfMissingFeatures(["STOKESIAN_DYNAMICS"]) +class StokesianDynamicsTest(sd.StokesianDynamicsTest): + device = 'cpu' + + def test_default(self): + self.falling_spheres(1.0, 1.0, 1.0) + + def test_rescaled(self): + self.falling_spheres(1.0, 4.5, 2.5) + + def test_different_time_step(self): + self.falling_spheres(0.7, 1.0, 1.0) + + def test_default_ft(self): + self.falling_spheres(1.0, 1.0, 1.0, 'ft') + + +@utx.skipIfMissingFeatures(["STOKESIAN_DYNAMICS"]) +class StokesianDiffusionTest(sd.StokesianDiffusionTest): + device = 'cpu' + + def test(self): + self.check() + + +if __name__ == '__main__': + ut.main() diff --git a/testsuite/python/stokesian_dynamics_gpu.py b/testsuite/python/stokesian_dynamics_gpu.py new file mode 100644 index 00000000000..b9ca6d0a73e --- /dev/null +++ b/testsuite/python/stokesian_dynamics_gpu.py @@ -0,0 +1,55 @@ +# +# Copyright (C) 2019-2020 The ESPResSo project +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +import unittest as ut +import unittest_decorators as utx +import stokesian_dynamics as sd + + +@utx.skipIfMissingGPU() +@utx.skipIfMissingFeatures(["STOKESIAN_DYNAMICS_GPU"]) +class StokesianDynamicsSetupTest(sd.StokesianDynamicsSetupTest): + device = 'gpu' + + def test_pbc_checks(self): + self.pbc_checks() + + +@utx.skipIfMissingGPU() +@utx.skipIfMissingFeatures(["STOKESIAN_DYNAMICS_GPU"]) +class StokesianDynamicsTest(sd.StokesianDynamicsTest): + device = 'gpu' + + def test_default_fts(self): + self.falling_spheres(1.0, 1.0, 1.0, 'fts', sd_short=True) + + def test_default_ft(self): + self.falling_spheres(1.0, 1.0, 1.0, 'ft', sd_short=True) + + +@utx.skipIfMissingGPU() +@utx.skipIfMissingFeatures(["STOKESIAN_DYNAMICS_GPU"]) +class StokesianDiffusionTest(sd.StokesianDiffusionTest): + device = 'gpu' + + def test(self): + pass + + +if __name__ == '__main__': + ut.main() diff --git a/testsuite/python/test_checkpoint.py b/testsuite/python/test_checkpoint.py index 4faf7b9d1f9..1ebba9fae27 100644 --- a/testsuite/python/test_checkpoint.py +++ b/testsuite/python/test_checkpoint.py @@ -35,6 +35,14 @@ and espressomd.has_features('ELECTROKINETICS')) +def skipIfMissingFeatureStokesianDynamics(): + """Specialized unittest skipIf decorator for missing Stokesian Dynamics.""" + if not espressomd.has_features(["STOKESIAN_DYNAMICS"]) and (not espressomd.has_features( + ["STOKESIAN_DYNAMICS_GPU"]) or not espressomd.gpu_available()): + return ut.skip("Skipping test: feature STOKESIAN_DYNAMICS unavailable") + return utx._id + + class CheckpointTest(ut.TestCase): @classmethod @@ -191,6 +199,14 @@ def test_thermostat_NPT(self): self.assertEqual(thmst['gamma0'], 2.0) self.assertEqual(thmst['gammav'], 0.1) + @skipIfMissingFeatureStokesianDynamics() + @ut.skipIf('THERM.SDM' not in modes, 'SDM thermostat not in modes') + def test_thermostat_SDM(self): + thmst = system.thermostat.get_state()[0] + self.assertEqual(thmst['type'], 'SD') + self.assertEqual(thmst['kT'], 1.0) + self.assertEqual(thmst['seed'], 42) + @utx.skipIfMissingFeatures('NPT') @ut.skipIf('INT.NPT' not in modes, 'NPT integrator not in modes') def test_integrator_NPT(self): @@ -233,6 +249,31 @@ def test_integrator_BD(self): params = integ.get_params() self.assertEqual(params, {}) + @utx.skipIfMissingFeatures('STOKESIAN_DYNAMICS') + @ut.skipIf('INT.SDM.CPU' not in modes, 'SDM CPU integrator not in modes') + def test_integrator_SDM_cpu(self): + integ = system.integrator.get_state() + self.assertIsInstance(integ, espressomd.integrate.StokesianDynamics) + expected_params = { + 'approximation_method': 'ft', 'device': 'cpu', 'radii': {0: 1.5}, + 'viscosity': 0.5, 'lubrication': False, 'pair_mobility': False, + 'self_mobility': True} + params = integ.get_params() + self.assertEqual(params, expected_params) + + @utx.skipIfMissingGPU() + @utx.skipIfMissingFeatures('STOKESIAN_DYNAMICS_GPU') + @ut.skipIf('INT.SDM.GPU' not in modes, 'SDM GPU integrator not in modes') + def test_integrator_SDM_gpu(self): + integ = system.integrator.get_state() + self.assertIsInstance(integ, espressomd.integrate.StokesianDynamics) + expected_params = { + 'approximation_method': 'fts', 'device': 'gpu', 'radii': {0: 1.0}, + 'viscosity': 2.0, 'lubrication': False, 'pair_mobility': True, + 'self_mobility': False} + params = integ.get_params() + self.assertEqual(params, expected_params) + @utx.skipIfMissingFeatures('LENNARD_JONES') @ut.skipIf('LJ' not in modes, "Skipping test due to missing mode.") def test_non_bonded_inter(self): diff --git a/testsuite/scripts/samples/CMakeLists.txt b/testsuite/scripts/samples/CMakeLists.txt index 69fcd4b459a..d9ca74c9e71 100644 --- a/testsuite/scripts/samples/CMakeLists.txt +++ b/testsuite/scripts/samples/CMakeLists.txt @@ -24,6 +24,10 @@ add_custom_target( sample_test(FILE test_billiard.py) sample_test(FILE test_chamber_game.py) sample_test(FILE test_constraints.py) +sample_test(FILE test_dancing.py SUFFIX cpu_ft) +sample_test(FILE test_dancing.py SUFFIX cpu_fts) +sample_test(FILE test_dancing.py SUFFIX gpu_ft LABELS "gpu") +sample_test(FILE test_dancing.py SUFFIX gpu_fts LABELS "gpu") sample_test(FILE test_diffusion_coefficient.py) sample_test(FILE test_dpd.py) sample_test(FILE test_drude_bmimpf6.py SUFFIX cpu) diff --git a/testsuite/scripts/samples/test_dancing.py b/testsuite/scripts/samples/test_dancing.py new file mode 100644 index 00000000000..1d322504438 --- /dev/null +++ b/testsuite/scripts/samples/test_dancing.py @@ -0,0 +1,63 @@ +# Copyright (C) 2020 The ESPResSo project +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +import unittest as ut +import importlib_wrapper +import numpy as np + +sd_device = "--gpu" if "gpu" in "@TEST_SUFFIX@" else "--cpu" +sd_method = "--fts" if "fts" in "@TEST_SUFFIX@" else "--ft" +if sd_method == "--fts" and sd_device == "--cpu": + y_min = -555 + intsteps = 5500 +else: + y_min = -200 + intsteps = 2100 + +sample, skipIfMissingFeatures = importlib_wrapper.configure_and_import( + "@SAMPLES_DIR@/dancing.py", script_suffix="@TEST_SUFFIX@", + gpu=(sd_device == "--gpu"), cmd_arguments=[sd_device, sd_method], + intsteps=intsteps, + ref_data="@CMAKE_SOURCE_DIR@/testsuite/python/data/dancing.txt") + + +@skipIfMissingFeatures +class Sample(ut.TestCase): + system = sample.system + + def test_trajectory(self): + # compare the simulated and published trajectories + simul = sample.positions[:, :, 0:2] + paper = sample.data.reshape([-1, 3, 2]) + + for pid in range(3): + dist = [] + # the simulated trajectory is oversampled by a ratio of 60:1 + # compared to the published trajectory (60 +/- 5 to 1) + for desired in paper[:, pid]: + if desired[1] < y_min: + break + # find the closest point in the simulated trajectory + idx = np.abs(simul[:, pid, 1] - desired[1]).argmin() + actual = simul[idx, pid] + dist.append(np.linalg.norm(actual - desired)) + self.assertLess(idx, sample.intsteps, msg='Insufficient sampling') + np.testing.assert_allclose(dist, 0, rtol=0, atol=0.7) + + +if __name__ == "__main__": + ut.main()