From f6a558f7f513d29dcbfbcb683a29800dbb17b599 Mon Sep 17 00:00:00 2001 From: "W. Kwabena Darko" Date: Tue, 3 Dec 2024 19:36:41 -0600 Subject: [PATCH 01/10] Implement MPCD reverse nonequilibrium shear flow on CPU --- hoomd/mpcd/CMakeLists.txt | 4 + hoomd/mpcd/ReverseNonequilibriumShearFlow.cc | 410 ++++++++++++++++++ hoomd/mpcd/ReverseNonequilibriumShearFlow.h | 119 +++++ .../ReverseNonequilibriumShearFlowUtilities.h | 119 +++++ hoomd/mpcd/__init__.py | 2 + hoomd/mpcd/module.cc | 2 + hoomd/mpcd/pytest/CMakeLists.txt | 1 + hoomd/mpcd/pytest/test_update.py | 207 +++++++++ hoomd/mpcd/update.py | 173 ++++++++ sphinx-doc/hoomd/module-mpcd.rst | 1 + sphinx-doc/hoomd/mpcd/module-update.rst | 13 + .../update/reversenonequilibriumshearflow.rst | 8 + sphinx-doc/module-mpcd-update.rst | 21 + 13 files changed, 1080 insertions(+) create mode 100644 hoomd/mpcd/ReverseNonequilibriumShearFlow.cc create mode 100644 hoomd/mpcd/ReverseNonequilibriumShearFlow.h create mode 100644 hoomd/mpcd/ReverseNonequilibriumShearFlowUtilities.h create mode 100644 hoomd/mpcd/pytest/test_update.py create mode 100644 hoomd/mpcd/update.py create mode 100644 sphinx-doc/hoomd/mpcd/module-update.rst create mode 100644 sphinx-doc/hoomd/mpcd/update/reversenonequilibriumshearflow.rst create mode 100644 sphinx-doc/module-mpcd-update.rst diff --git a/hoomd/mpcd/CMakeLists.txt b/hoomd/mpcd/CMakeLists.txt index 40b5ce9ac8..496a125a95 100644 --- a/hoomd/mpcd/CMakeLists.txt +++ b/hoomd/mpcd/CMakeLists.txt @@ -15,6 +15,7 @@ set(_mpcd_cc_sources ManualVirtualParticleFiller.cc ParallelPlateGeometryFiller.cc PlanarPoreGeometryFiller.cc + ReverseNonequilibriumShearFlow.cc Sorter.cc SRDCollisionMethod.cc StreamingMethod.cc @@ -41,6 +42,8 @@ set(_mpcd_headers ParallelPlateGeometryFiller.h PlanarPoreGeometryFiller.h RejectionVirtualParticleFiller.h + ReverseNonequilibriumShearFlow.h + ReverseNonequilibriumShearFlowUtilities.h Sorter.h SRDCollisionMethod.h StreamingMethod.h @@ -265,6 +268,7 @@ set(files methods.py stream.py tune.py + update.py ) install(FILES ${files} DESTINATION ${PYTHON_SITE_INSTALL_DIR}/mpcd) copy_files_to_build("${files}" "mpcd" "*.py") diff --git a/hoomd/mpcd/ReverseNonequilibriumShearFlow.cc b/hoomd/mpcd/ReverseNonequilibriumShearFlow.cc new file mode 100644 index 0000000000..05130381d0 --- /dev/null +++ b/hoomd/mpcd/ReverseNonequilibriumShearFlow.cc @@ -0,0 +1,410 @@ +// Copyright (c) 2009-2024 The Regents of the University of Michigan. +// Part of HOOMD-blue, released under the BSD 3-Clause License. + +/*! + * \file ReverseNonequilibriumShearFlow.cc + * \brief Definition of reverse nonequilibrium shear flow updater + */ + +#include + +#include "ReverseNonequilibriumShearFlow.h" +#include "ReverseNonequilibriumShearFlowUtilities.h" + +namespace hoomd + { +/*! + * \param sysdata MPCD system data this updater will act on + * \param num_swap Max number of swaps + * \param slab_width Slab width + * \param target_momentum Target momentum in x-direction for the two slabs + */ +mpcd::ReverseNonequilibriumShearFlow::ReverseNonequilibriumShearFlow( + std::shared_ptr sysdef, + std::shared_ptr trigger, + unsigned int num_swap, + Scalar slab_width, + Scalar target_momentum) + : Updater(sysdef, trigger), m_mpcd_pdata(sysdef->getMPCDParticleData()), m_num_swap(num_swap), + m_slab_width(slab_width), m_summed_momentum_exchange(0), m_num_lo(0), + m_particles_lo(m_exec_conf), m_num_hi(0), m_particles_hi(m_exec_conf), + m_target_momentum(target_momentum), m_update_slabs(true) + { + m_exec_conf->msg->notice(5) << "Constructing ReverseNonequilibriumShearFlow" << std::endl; + + m_pdata->getBoxChangeSignal() + .connect(this); + + GPUArray particles_staged(2 * m_num_swap, m_exec_conf); + m_particles_staged.swap(particles_staged); + } + +mpcd::ReverseNonequilibriumShearFlow::~ReverseNonequilibriumShearFlow() + { + m_exec_conf->msg->notice(5) << "Destroying ReverseNonequilibriumShearFlow" << std::endl; + m_pdata->getBoxChangeSignal() + .disconnect(this); + } + +/*! + * \param num_swap Max number of swaps + */ +void mpcd::ReverseNonequilibriumShearFlow::setNumSwap(unsigned int num_swap) + { + if (num_swap > m_num_swap) + { + GPUArray particles_staged(2 * num_swap, m_exec_conf); + m_particles_staged.swap(particles_staged); + } + m_num_swap = num_swap; + } + +//! Set the target momentum +void mpcd::ReverseNonequilibriumShearFlow::setTargetMomentum(Scalar target_momentum) + { + m_target_momentum = target_momentum; + } + +/*! + * \param slab_width Slab width + */ +void mpcd::ReverseNonequilibriumShearFlow::setSlabWidth(Scalar slab_width) + { + m_slab_width = slab_width; + requestUpdateSlabs(); + } + +void mpcd::ReverseNonequilibriumShearFlow::setSlabs() + { + const BoxDim& global_box = m_pdata->getGlobalBox(); + if (m_slab_width > Scalar(0.5) * global_box.getL().y) + { + throw std::runtime_error("Slab width cannot be larger than Ly/2"); + } + + const Scalar3 global_lo = global_box.getLo(); + m_pos_lo = make_scalar2(global_lo.y, global_lo.y + m_slab_width); + m_pos_hi = make_scalar2(Scalar(0.0), m_slab_width); + } + +/*! + * \param timestep Current time step of the simulation + */ +void mpcd::ReverseNonequilibriumShearFlow::update(uint64_t timestep) + { + // if slabs have changed, update them + if (m_update_slabs) + { + setSlabs(); + m_update_slabs = false; + } + + findSwapParticles(); + stageSwapParticles(); + swapParticleMomentum(); + } + +/*! + * Finds all particles in the "lower" and "upper" slab in y direction and + * puts them into two GPUArrays, sorted by their momentum closest to -/+ target_momentum in x + * direction. + */ +void mpcd::ReverseNonequilibriumShearFlow::findSwapParticles() + { + // fill the layers with (unsorted) particles. this uses do loop for reallocation + bool filled = false; + do + { + const size_t num_lo_alloc = m_particles_lo.getNumElements(); + const size_t num_hi_alloc = m_particles_hi.getNumElements(); + { + ArrayHandle h_pos(m_mpcd_pdata->getPositions(), + access_location::host, + access_mode::read); + ArrayHandle h_vel(m_mpcd_pdata->getVelocities(), + access_location::host, + access_mode::read); + + ArrayHandle h_particles_lo(m_particles_lo, + access_location::host, + access_mode::overwrite); + ArrayHandle h_particles_hi(m_particles_hi, + access_location::host, + access_mode::overwrite); + + // filter particles into their slab in y-direction and record momentum in x-direction + m_num_lo = 0; + m_num_hi = 0; + for (unsigned int idx = 0; idx < m_mpcd_pdata->getN(); ++idx) + { + const Scalar4 vel = h_vel.data[idx]; + const Scalar momentum = vel.x * m_mpcd_pdata->getMass(); + const Scalar y = h_pos.data[idx].y; + if (m_pos_lo.x <= y && y < m_pos_lo.y + && momentum > Scalar(0.0)) // lower slab, search for positive momentum + { + if (m_num_lo < num_lo_alloc) + { + h_particles_lo.data[m_num_lo] + = make_scalar2(__int_as_scalar(idx), momentum); + } + ++m_num_lo; + } + else if (m_pos_hi.x <= y && y < m_pos_hi.y + && momentum < Scalar(0.0)) // higher slab, search for negative momentum + { + if (m_num_hi < num_hi_alloc) + { + h_particles_hi.data[m_num_hi] + = make_scalar2(__int_as_scalar(idx), momentum); + } + ++m_num_hi; + } + } + } + filled = true; + + // reallocate if required and go again, this won't happen too much + if (m_num_lo > num_lo_alloc) + { + GPUArray particles_lo(m_num_lo, m_exec_conf); + m_particles_lo.swap(particles_lo); + filled = false; + } + if (m_num_hi > num_hi_alloc) + { + GPUArray particles_hi(m_num_hi, m_exec_conf); + m_particles_hi.swap(particles_hi); + filled = false; + } + + } while (!filled); + + // sort the arrays + { + ArrayHandle h_particles_lo(m_particles_lo, + access_location::host, + access_mode::readwrite); + ArrayHandle h_particles_hi(m_particles_hi, + access_location::host, + access_mode::readwrite); + ArrayHandle h_tag(m_mpcd_pdata->getTags(), + access_location::host, + access_mode::read); + + if (std::isinf(m_target_momentum)) + { + std::sort(h_particles_lo.data, + h_particles_lo.data + m_num_lo, + detail::MaximumMomentum(h_tag.data)); + std::sort(h_particles_hi.data, + h_particles_hi.data + m_num_hi, + detail::MinimumMomentum(h_tag.data)); + } + else + { + std::sort(h_particles_lo.data, + h_particles_lo.data + m_num_lo, + detail::CompareMomentumToTarget(m_target_momentum, h_tag.data)); + std::sort(h_particles_hi.data, + h_particles_hi.data + m_num_hi, + detail::CompareMomentumToTarget(-m_target_momentum, h_tag.data)); + } + } + } + +/*! + * Stage particles momenta from both slabs into a queue for swapping + */ +void mpcd::ReverseNonequilibriumShearFlow::stageSwapParticles() + { + // determine number of pairs to swap + unsigned int num_lo_global = m_num_lo; + unsigned int num_hi_global = m_num_hi; +#ifdef ENABLE_MPI + if (m_sysdef->isDomainDecomposed()) + { + auto mpi_comm = m_exec_conf->getMPICommunicator(); + MPI_Allreduce(MPI_IN_PLACE, &num_lo_global, 1, MPI_UNSIGNED, MPI_SUM, mpi_comm); + MPI_Allreduce(MPI_IN_PLACE, &num_hi_global, 1, MPI_UNSIGNED, MPI_SUM, mpi_comm); + } +#endif // ENABLE_MPI + const unsigned int num_pairs = std::min(m_num_swap, std::min(num_lo_global, num_hi_global)); + + // stage swaps into queue + ArrayHandle h_particles_lo(m_particles_lo, access_location::host, access_mode::read); + ArrayHandle h_particles_hi(m_particles_hi, access_location::host, access_mode::read); + ArrayHandle h_particles_staged(m_particles_staged, + access_location::host, + access_mode::overwrite); + m_num_staged = 0; + unsigned int lo_idx = 0; + unsigned int hi_idx = 0; + for (unsigned int i = 0; i < num_pairs; ++i) + { +#ifdef ENABLE_MPI + if (m_sysdef->isDomainDecomposed()) + { + const auto mpi_comm = m_exec_conf->getMPICommunicator(); + const int rank = m_exec_conf->getRank(); + + // find particle in lo slab from all ranks (p > 0) + MPI_Op lo_op; + Scalar_Int lo_swap; + lo_swap.i = rank; + if (std::isinf(m_target_momentum)) + { + lo_op = MPI_MAXLOC; + lo_swap.s = (lo_idx < m_num_lo) ? h_particles_lo.data[lo_idx].y : Scalar(0.0); + } + else + { + lo_op = MPI_MINLOC; + lo_swap.s = (lo_idx < m_num_lo) + ? std::fabs(h_particles_lo.data[lo_idx].y - m_target_momentum) + : std::numeric_limits::max(); + } + MPI_Allreduce(MPI_IN_PLACE, &lo_swap, 1, MPI_HOOMD_SCALAR_INT, lo_op, mpi_comm); + + // find particle in hi slab to swap (p < 0) + Scalar_Int hi_swap; + hi_swap.i = rank; + if (hi_idx < m_num_hi) + { + hi_swap.s = (std::isinf(m_target_momentum)) + ? h_particles_hi.data[hi_idx].y + : std::fabs(h_particles_hi.data[hi_idx].y + m_target_momentum); + } + else + { + hi_swap.s = std::numeric_limits::max(); + } + MPI_Allreduce(MPI_IN_PLACE, &hi_swap, 1, MPI_HOOMD_SCALAR_INT, MPI_MINLOC, mpi_comm); + + if (lo_swap.i == rank || hi_swap.i == rank) // at most, 2 ranks participate in the swap + { + if (lo_swap.i != hi_swap.i) // particles to swap are on different ranks, needs MPI + { + Scalar2 particle; + unsigned int dest; + if (lo_swap.i == rank) + { + particle = h_particles_lo.data[lo_idx++]; + dest = hi_swap.i; + } + else + { + particle = h_particles_hi.data[hi_idx++]; + dest = lo_swap.i; + } + + MPI_Sendrecv_replace(&particle.y, + 1, + MPI_HOOMD_SCALAR, + dest, + 0, + dest, + 0, + mpi_comm, + MPI_STATUS_IGNORE); + + h_particles_staged.data[m_num_staged++] = particle; + } + else // particles are on the same rank + { + Scalar2 lo_particle = h_particles_lo.data[lo_idx++]; + Scalar2 hi_particle = h_particles_hi.data[hi_idx++]; + std::swap(lo_particle.y, hi_particle.y); + h_particles_staged.data[m_num_staged++] = lo_particle; + h_particles_staged.data[m_num_staged++] = hi_particle; + } + } + } + else +#endif // ENABLE_MPI + { + Scalar2 lo_particle = h_particles_lo.data[lo_idx++]; + Scalar2 hi_particle = h_particles_hi.data[hi_idx++]; + std::swap(lo_particle.y, hi_particle.y); + h_particles_staged.data[m_num_staged++] = lo_particle; + h_particles_staged.data[m_num_staged++] = hi_particle; + } + } + } + +/*! + * Apply new momenta to particles from the queue. + */ +void mpcd::ReverseNonequilibriumShearFlow::swapParticleMomentum() + { + ArrayHandle h_vel(m_mpcd_pdata->getVelocities(), + access_location::host, + access_mode::readwrite); + ArrayHandle h_particles_staged(m_particles_staged, + access_location::host, + access_mode::read); + + // perform swap and sum momentum exchange + Scalar momentum_sum(0); + const Scalar mass = m_mpcd_pdata->getMass(); + for (unsigned int i = 0; i < m_num_staged; ++i) + { + const Scalar2 pidx_mom = h_particles_staged.data[i]; + const unsigned int pidx = __scalar_as_int(pidx_mom.x); + const Scalar new_momentum = pidx_mom.y; + + const Scalar current_momentum = h_vel.data[pidx].x * mass; + + h_vel.data[pidx].x = new_momentum / mass; + momentum_sum += std::fabs(new_momentum - current_momentum); + } + +#ifdef ENABLE_MPI + if (m_sysdef->isDomainDecomposed()) + { + auto comm = m_exec_conf->getMPICommunicator(); + MPI_Allreduce(MPI_IN_PLACE, &momentum_sum, 1, MPI_HOOMD_SCALAR, MPI_SUM, comm); + } +#endif + + // absolute value gives extra factor of 2, divide out when accumulating + m_summed_momentum_exchange += Scalar(0.5) * momentum_sum; + } + +namespace mpcd + { +namespace detail + { +/*! + * \param m Python module to export to + */ +void export_ReverseNonequilibriumShearFlow(pybind11::module& m) + { + namespace py = pybind11; + py::class_>( + m, + "ReverseNonequilibriumShearFlow") + .def(py::init, + std::shared_ptr, + unsigned int, + Scalar, + Scalar>()) + .def_property("num_swaps", + &mpcd::ReverseNonequilibriumShearFlow::getNumSwap, + &mpcd::ReverseNonequilibriumShearFlow::setNumSwap) + .def_property("slab_width", + &mpcd::ReverseNonequilibriumShearFlow::getSlabWidth, + &mpcd::ReverseNonequilibriumShearFlow::setSlabWidth) + .def_property("target_momentum", + &mpcd::ReverseNonequilibriumShearFlow::getTargetMomentum, + &mpcd::ReverseNonequilibriumShearFlow::setTargetMomentum) + .def_property_readonly("summed_exchanged_momentum", + &mpcd::ReverseNonequilibriumShearFlow::getSummedExchangedMomentum); + } + } // namespace detail + } // namespace mpcd + } // end namespace hoomd diff --git a/hoomd/mpcd/ReverseNonequilibriumShearFlow.h b/hoomd/mpcd/ReverseNonequilibriumShearFlow.h new file mode 100644 index 0000000000..c56743438a --- /dev/null +++ b/hoomd/mpcd/ReverseNonequilibriumShearFlow.h @@ -0,0 +1,119 @@ +// Copyright (c) 2009-2024 The Regents of the University of Michigan. +// Part of HOOMD-blue, released under the BSD 3-Clause License. + +/*! + * \file ReverseNonequilibriumShearFlow.h + * \brief Declaration of Reverse nonequilibrium shear flow + */ + +#ifndef MPCD_REVERSE_NONEQUILIBRIUM_SHEAR_FLOW_H_ +#define MPCD_REVERSE_NONEQUILIBRIUM_SHEAR_FLOW_H_ + +#ifdef NVCC +#error This header cannot be compiled by nvcc +#endif + +#include "hoomd/ParticleGroup.h" +#include "hoomd/Updater.h" +#include + +namespace hoomd + { +namespace mpcd + { +//! Reverse nonequilibrium shear flow updater +/*! + * A flow is induced by swapping velocities in x direction based on particle position in + * y-direction. + */ +class PYBIND11_EXPORT ReverseNonequilibriumShearFlow : public Updater + { + public: + //! Constructor + ReverseNonequilibriumShearFlow(std::shared_ptr sysdef, + std::shared_ptr trigger, + unsigned int num_swap, + Scalar slab_width, + Scalar target_momentum); + + //! Destructor + virtual ~ReverseNonequilibriumShearFlow(); + + //! Apply velocity swaps + virtual void update(uint64_t timestep); + + //! Get max number of swaps + Scalar getNumSwap() const + { + return m_num_swap; + } + + //! Set the maximum number of swapped pairs + void setNumSwap(unsigned int num_swap); + + //! Get slab width + Scalar getSlabWidth() const + { + return m_slab_width; + } + + //! Set the slab width + void setSlabWidth(Scalar slab_width); + + //! Get target momentum + Scalar getTargetMomentum() const + { + return m_target_momentum; + } + + //! Set the target momentum + void setTargetMomentum(Scalar target_momentum); + + //! Get summed exchanged momentum + Scalar getSummedExchangedMomentum() const + { + return m_summed_momentum_exchange; + } + + protected: + std::shared_ptr m_mpcd_pdata; //!< MPCD particle data + unsigned int m_num_swap; //!< maximum number of swaps + Scalar m_slab_width; //!< width of slabs + Scalar m_summed_momentum_exchange; //!< summed momentum excange between slabs + Scalar2 m_pos_lo; //!< position of bottom slab in box + Scalar2 m_pos_hi; //!< position of top slab in box + unsigned int m_num_lo; //!< number of particles in bottom slab + GPUArray m_particles_lo; //!< List of all particles (indices,momentum) in bottom slab + //!< sorted by momentum closest to +m_target_momentum + unsigned int m_num_hi; //!< number of particles in top slab + GPUArray m_particles_hi; //!< List of all particles (indices,momentum) in top slab + //!< sorted by momentum closest to -m_target_momentum + unsigned int m_num_staged; //!< number of particles staged for swapping + GPUArray m_particles_staged; //!< List of all particles staged for swapping + Scalar m_target_momentum; //!< target momentum for particles in the slabs + + //! Find candidate particles for swapping in the slabs + virtual void findSwapParticles(); + + //! Stage particle momentum for swapping + void stageSwapParticles(); + + //! Swaps momentum between the slabs + virtual void swapParticleMomentum(); + + private: + bool m_update_slabs; //!< If true, update the slab positions + + //! Request to check box on next update + void requestUpdateSlabs() + { + m_update_slabs = true; + } + + //! Sets the slab positions in the box and validates them + void setSlabs(); + }; + + } // end namespace mpcd + } // end namespace hoomd +#endif // MPCD_REVERSE_NONEQUILIBRIUM_SHEAR_FLOW_H_ diff --git a/hoomd/mpcd/ReverseNonequilibriumShearFlowUtilities.h b/hoomd/mpcd/ReverseNonequilibriumShearFlowUtilities.h new file mode 100644 index 0000000000..f88a8184cb --- /dev/null +++ b/hoomd/mpcd/ReverseNonequilibriumShearFlowUtilities.h @@ -0,0 +1,119 @@ +// Copyright (c) 2009-2024 The Regents of the University of Michigan. +// Part of HOOMD-blue, released under the BSD 3-Clause License. + +/*! + * \file ReverseNonequilibriumShearFlowUtilities.h + * \brief Helper functions for Reverse nonequilibrium shear flow + */ + +#ifndef MPCD_REVERSE_NONEQUILIBRIUM_SHEAR_FLOW_UTILITIES_H_ +#define MPCD_REVERSE_NONEQUILIBRIUM_SHEAR_FLOW_UTILITIES_H_ + +#include "hoomd/HOOMDMath.h" + +namespace hoomd + { +namespace mpcd + { +namespace detail + { + +class CompareMomentumToTarget + { + public: + CompareMomentumToTarget(Scalar target_momentum, const unsigned int* tags_) + : p(target_momentum), tags(tags_) + { + } + + bool operator()(const Scalar2& in0, const Scalar2& in1) const + { + const Scalar dp0 = std::fabs(in0.y - p); + const Scalar dp1 = std::fabs(in1.y - p); + if (dp0 < dp1) + { + // 0 is closer to target than 1 + return true; + } + else if (dp0 > dp1) + { + // 0 is farther from target than 1 + return false; + } + else + { + // both are equal distant, break tie using lower tag + return tags[__scalar_as_int(in0.x)] < tags[__scalar_as_int(in1.x)]; + } + } + + private: + const Scalar p; //!< Momentum target + const unsigned int* const tags; //!< Tag array + }; + +class MaximumMomentum + { + public: + MaximumMomentum(const unsigned int* tags_) : tags(tags_) { } + + bool operator()(const Scalar2& in0, const Scalar2& in1) const + { + const Scalar p0 = in0.y; + const Scalar p1 = in1.y; + if (p0 > p1) + { + // particle 0 has higher momemtum than 1 so should be selected first + return true; + } + else if (p0 < p1) + { + // particle 0 has lower momentum than 1 + return false; + } + else + { + // both are equal distant, break tie using lower tag + return tags[__scalar_as_int(in0.x)] < tags[__scalar_as_int(in1.x)]; + } + } + + private: + const unsigned int* const tags; //!< Tag array + }; + +class MinimumMomentum + { + public: + MinimumMomentum(const unsigned int* tags_) : tags(tags_) { } + + bool operator()(const Scalar2& in0, const Scalar2& in1) const + { + const Scalar p0 = in0.y; + const Scalar p1 = in1.y; + if (p0 < p1) + { + // particle 0 is lower in momemtum than 1 so should be selected first + return true; + } + else if (p0 > p1) + { + // partilce 0 has higher momentum than 1 + return false; + } + else + { + // both are equal distant, break tie using lower tag + return tags[__scalar_as_int(in0.x)] < tags[__scalar_as_int(in1.x)]; + } + } + + private: + const unsigned int* const tags; //!< Tag array + }; + + } // end namespace detail + } // end namespace mpcd + } // end namespace hoomd + +#endif // MPCD_REVERSE_NONEQUILIBRIUM_SHEAR_FLOW_UTILITIES_H_ diff --git a/hoomd/mpcd/__init__.py b/hoomd/mpcd/__init__.py index fe9eae71f8..e32dc11758 100644 --- a/hoomd/mpcd/__init__.py +++ b/hoomd/mpcd/__init__.py @@ -67,6 +67,7 @@ from hoomd.mpcd import methods from hoomd.mpcd import stream from hoomd.mpcd import tune +from hoomd.mpcd import update __all__ = [ "Integrator", @@ -77,4 +78,5 @@ "methods", "stream", "tune", + "update", ] diff --git a/hoomd/mpcd/module.cc b/hoomd/mpcd/module.cc index 92b71cf59e..d97d8bb400 100644 --- a/hoomd/mpcd/module.cc +++ b/hoomd/mpcd/module.cc @@ -42,6 +42,7 @@ void export_ParallelPlateGeometry(pybind11::module&); void export_ParallelPlateGeometryFiller(pybind11::module&); void export_PlanarPoreGeometry(pybind11::module&); void export_PlanarPoreGeometryFiller(pybind11::module&); +void export_ReverseNonequilibriumShearFlow(pybind11::module&); void export_Sorter(pybind11::module&); void export_SphereGeometry(pybind11::module&); void export_SphereGeometryFiller(pybind11::module&); @@ -218,6 +219,7 @@ PYBIND11_MODULE(_mpcd, m) export_ParallelPlateGeometryFiller(m); export_PlanarPoreGeometry(m); export_PlanarPoreGeometryFiller(m); + export_ReverseNonequilibriumShearFlow(m); export_Sorter(m); export_SphereGeometry(m); export_SphereGeometryFiller(m); diff --git a/hoomd/mpcd/pytest/CMakeLists.txt b/hoomd/mpcd/pytest/CMakeLists.txt index cbdddfcaae..400dd0860a 100644 --- a/hoomd/mpcd/pytest/CMakeLists.txt +++ b/hoomd/mpcd/pytest/CMakeLists.txt @@ -9,6 +9,7 @@ set(files __init__.py test_snapshot.py test_stream.py test_tune.py + test_update.py ) install(FILES ${files} diff --git a/hoomd/mpcd/pytest/test_update.py b/hoomd/mpcd/pytest/test_update.py new file mode 100644 index 0000000000..5919ac0c74 --- /dev/null +++ b/hoomd/mpcd/pytest/test_update.py @@ -0,0 +1,207 @@ +# Copyright (c) 2009-2024 The Regents of the University of Michigan. +# Part of HOOMD-blue, released under the BSD 3-Clause License. + +import pytest +import hoomd +import numpy as np + +from hoomd.conftest import pickling_check +from hoomd.logging import LoggerCategories +from hoomd.conftest import logging_check + + +@pytest.fixture +def snap(): + snap = hoomd.Snapshot() + if snap.communicator.rank == 0: + snap.configuration.box = [20, 20, 20, 0, 0, 0] + snap.particles.N = 0 + snap.particles.types = ["A"] + snap.mpcd.types = ["A"] + snap.mpcd.mass = 1.0 + return snap + + +@pytest.fixture +def flow(): + flow = hoomd.mpcd.update.ReverseNonequilibriumShearFlow( + trigger=1, num_swaps=1, slab_width=1.0 + ) + return flow + + +decompositions = [None] +if hoomd.communicator.Communicator().num_ranks == 2: + decompositions += [(1, 2, 1), (2, 1, 1)] + + +class TestReverseNonequilibriumShearFlow: + def test_create(self, simulation_factory, snap, flow): + # before attachment + flow.num_swaps = 5 # changing the value before attachment + flow.target_momentum = 2.5 # changing the value before attachment + trigger = hoomd.trigger.Periodic(1) + flow.trigger = trigger + + assert flow.trigger is trigger + assert flow.slab_width == 1.0 + + with pytest.raises(hoomd.error.DataAccessError): + assert flow.summed_exchanged_momentum == 0.0 + + sim = simulation_factory(snap) + sim.operations.integrator = hoomd.mpcd.Integrator(dt=0.09) + sim.operations.updaters.append(flow) + sim.run(0) + + # after attachment + assert flow.num_swaps == 5 + assert flow.target_momentum == 2.5 + assert flow.trigger is trigger + assert flow.slab_width == 1.0 + + assert flow.summed_exchanged_momentum == 0.0 + + def test_pickling(self, simulation_factory, snap, flow): + pickling_check(flow) + + flow.trigger = hoomd.trigger.Periodic(1) + flow.num_swaps = 5 + pickling_check(flow) + + sim = simulation_factory(snap) + sim.operations.integrator = hoomd.mpcd.Integrator(dt=0.09) + sim.operations.updaters.append(flow) + sim.run(0) + pickling_check(flow) + + # particle locations in/out of the slabs; and at the boundaries + @pytest.mark.parametrize("decomposition", decompositions) + @pytest.mark.parametrize( + "v0, v1, y0, y1, swap", + [ + [-2, 2, 0.5, -9.5, True], + [2, -2, 0.5, -9.5, False], + [-2, -2, 0.5, -9.5, False], + [2, 2, 0.5, -9.5, False], + [-2, 2, 0.5, -7.0, False], + [-2, 2, 5.0, -9.5, False], + [-2, 2, 5.0, -7.0, False], + [-2, 2, 0.0, -10.0, True], + [-2, 2, 0.0, -9.0, False], + [-2, 2, 1.0, -9.0, False], + [-2, 2, 1.0, -10.0, False], + ], + ) + def test_particle_swapping( + self, simulation_factory, snap, flow, v0, v1, y0, y1, swap, decomposition + ): + if snap.communicator.rank == 0: + snap.mpcd.N = 2 + snap.mpcd.position[0] = [1.0, y0, 2.0] + snap.mpcd.position[1] = [1.0, y1, 2.0] + snap.mpcd.velocity[0] = [v0, 4.0, 5.0] + snap.mpcd.velocity[1] = [v1, 1.0, 3.0] + + flow.target_momentum = 2.0 + sim = simulation_factory(snap, domain_decomposition=decomposition) + sim.operations.integrator = hoomd.mpcd.Integrator(dt=0.09) + sim.operations.updaters.append(flow) + sim.run(1) + + snap = sim.state.get_snapshot() + if snap.communicator.rank == 0: + assert np.allclose(snap.mpcd.velocity[0], [v1 if swap else v0, 4.0, 5.0]) + assert np.allclose(snap.mpcd.velocity[1], [v0 if swap else v1, 1.0, 3.0]) + + sim.run(1) + + # confirm that no further swaps happen after the previous run + snapshot = sim.state.get_snapshot() + if snapshot.communicator.rank == 0: + assert np.allclose(snapshot.mpcd.velocity[0], snap.mpcd.velocity[0]) + assert np.allclose(snapshot.mpcd.velocity[1], snap.mpcd.velocity[1]) + + def test_swapping_at_different_timesteps(self, simulation_factory, snap, flow): + flow.target_momentum = 2.0 + # set initial positions and velocities + if snap.communicator.rank == 0: + snap.mpcd.N = 4 + snap.mpcd.position[0] = [1.0, 0.0, 2.0] + snap.mpcd.position[1] = [2.0, 0.3, 1.0] + snap.mpcd.position[2] = [1.0, -9.7, 3.0] + snap.mpcd.position[3] = [3.0, -10.0, 1.0] + + snap.mpcd.velocity[0] = [-1.0, 4.0, 5.0] + snap.mpcd.velocity[1] = [-9.0, 1.0, 3.0] + snap.mpcd.velocity[2] = [1.0, -2.0, 4.0] + snap.mpcd.velocity[3] = [5.0, -4.0, -1.0] + + # set up simulation + flow.num_swaps = 1 + sim = simulation_factory(snap) + sim.operations.integrator = hoomd.mpcd.Integrator(dt=0.09) + sim.operations.updaters.append(flow) + sim.run(1) + + assert flow.summed_exchanged_momentum == 2.0 + + snap = sim.state.get_snapshot() + if snap.communicator.rank == 0: + # check that swap only happens between particle ids 0,2 + assert np.allclose(snap.mpcd.velocity[0], [1.0, 4.0, 5.0]) + assert np.allclose(snap.mpcd.velocity[1], [-9.0, 1.0, 3.0]) + assert np.allclose(snap.mpcd.velocity[2], [-1.0, -2.0, 4.0]) + assert np.allclose(snap.mpcd.velocity[3], [5.0, -4.0, -1.0]) + + sim.run(1) + + assert flow.summed_exchanged_momentum == 2.0 + 14.0 + + snap = sim.state.get_snapshot() + if snap.communicator.rank == 0: + # check that swapping is now between particle ids 1,3 + assert np.allclose(snap.mpcd.velocity[0], [1.0, 4.0, 5.0]) + assert np.allclose(snap.mpcd.velocity[1], [5.0, 1.0, 3.0]) + assert np.allclose(snap.mpcd.velocity[2], [-1.0, -2.0, 4.0]) + assert np.allclose(snap.mpcd.velocity[3], [-9.0, -4.0, -1.0]) + + def test_default_target_momentum(self, simulation_factory, snap, flow): + if snap.communicator.rank == 0: + snap.mpcd.N = 4 + # set initial positions and velocities + snap.mpcd.position[0] = [1.0, 0.5, 2.0] + snap.mpcd.position[1] = [2.0, 0.3, 1.0] + snap.mpcd.position[2] = [1.0, -9.7, 3.0] + snap.mpcd.position[3] = [3.0, -9.2, 1.0] + + snap.mpcd.velocity[0] = [-9.0, 4.0, 5.0] + snap.mpcd.velocity[1] = [-1.0, 1.0, 3.0] + snap.mpcd.velocity[2] = [1.0, -2.0, 4.0] + snap.mpcd.velocity[3] = [3.0, -4.0, -1.0] + + # set up simulation + sim = simulation_factory(snap) + sim.operations.integrator = hoomd.mpcd.Integrator(dt=0.09) + sim.operations.updaters.append(flow) + sim.run(1) + + snap = sim.state.get_snapshot() + if snap.communicator.rank == 0: + # check that swap only happens between fastest and slowest particles + assert np.allclose(snap.mpcd.velocity[0], [3.0, 4.0, 5.0]) + assert np.allclose(snap.mpcd.velocity[1], [-1.0, 1.0, 3.0]) + assert np.allclose(snap.mpcd.velocity[2], [1.0, -2.0, 4.0]) + assert np.allclose(snap.mpcd.velocity[3], [-9.0, -4.0, -1.0]) + + def test_logging(self): + logging_check( + hoomd.mpcd.update.ReverseNonequilibriumShearFlow, + ("mpcd", "update"), + { + "summed_exchanged_momentum": { + "category": LoggerCategories.scalar, + "default": True, + } + }, + ) diff --git a/hoomd/mpcd/update.py b/hoomd/mpcd/update.py new file mode 100644 index 0000000000..296866ee40 --- /dev/null +++ b/hoomd/mpcd/update.py @@ -0,0 +1,173 @@ +# Copyright (c) 2009-2024 The Regents of the University of Michigan. +# Part of HOOMD-blue, released under the BSD 3-Clause License. + +r"""MPCD updaters. + +.. invisible-code-block: python + + simulation = hoomd.util.make_example_simulation(mpcd_types=["A"]) + simulation.operations.integrator = hoomd.mpcd.Integrator(dt=0.1) + +""" + +from . import _mpcd +from hoomd.operation import Updater +from hoomd.data.parameterdicts import ParameterDict +from hoomd.logging import log +from hoomd.data.typeconverter import OnlyTypes, positive_real + +import math + + +class ReverseNonequilibriumShearFlow(Updater): + r"""MPCD reverse nonequilibrium shear flow. + + Args: + trigger (hoomd.trigger.trigger_like): Trigger to swap momentum. + + num_swaps (int): Maximum number of times to swap momentum per update. + + slab_width (float): Width of momentum-exchange slabs. + + target_momentum (float): Target momentum for swapped particles. This + argument has a default value of infinity but can be + redefined with positive real numbers only. + + This updater generates a bidirectional shear flow in *x* by imposing a + momentum flux on the system in *y*. + There are two exchange slabs separated by a distance of :math:`L_y/2`. The + edges of these slabs are located at (:math:`-L_y/2`, :math:`-L_y/2` + + `slab_width`) and (:math:`0.0`, `slab_width`) along the *y*-direction + (gradient direction) of the simulation box. On each `trigger`, particles + are sorted into the exchange slabs based on their positions in the box. + Particles whose *x*-component momenta are near `target_momentum` + in the lower slab, and those near -`target_momentum` in the upper + slab, are selected for a pairwise momentum swap. Up to `num_swaps` swaps + are executed per update. + + The amount of momentum transferred from the lower slab to the upper slab is + known. Therefore `summed_exchanged_momentum`, which returns the accumulated + momentum exchanged till the current timestep, is used to calculate the + momentum flux, :math:`j_{xy}`. The shear rate, :math:`\dot{\gamma}`, is + also extracted as the gradient of the linear velocity profile developed + from the flow. The viscosity can then be computed from these two quantities + as: + + .. math:: + + \eta (\dot{\gamma}) = \frac{j_{xy}}{\dot{\gamma}}. + + .. rubric:: Examples: + + In the original implementation by + `Müller-Plathe `_, + only the fastest particle and the slowest particle are swapped. This is + achieved by setting `num_swaps` to 1 and keeping `target_momentum` at its + default value of infinity. + + .. code-block:: python + + flow = hoomd.mpcd.update.ReverseNonequilibriumShearFlow( + trigger=1, num_swaps=1, slab_width=1 + ) + simulation.operations.updaters.append(flow) + + An alternative approach proposed by + `Tenney and Maginn `_ swaps particles + that are instead closest to the `target_momentum`, typically requiring more + swaps per update. + + .. code-block:: python + + flow = hoomd.mpcd.update.ReverseNonequilibriumShearFlow( + trigger=1, num_swaps=10, slab_width=1, target_momentum=5 + ) + simulation.operations.updaters.append(flow) + + {inherited} + + ---------- + + **Members defined in** `ReverseNonequilibriumShearFlow`: + + Attributes: + trigger (hoomd.trigger.trigger_like): Trigger to swap momentum. + + .. rubric:: Example: + + .. code-block:: python + + flow.trigger = 1 + + num_swaps (int): Maximum number of times to swap momentum per update. + + .. rubric:: Example: + + .. code-block:: python + + flow.num_swaps = 10 + + slab_width (float): Width of momentum-exchange slabs. + + .. rubric:: Example: + + .. code-block:: python + + flow.slab_width = 1 + + target_momentum (float): Target momentum for swapped particles. + + .. rubric:: Example: + + .. code-block:: python + + flow.target_momentum = 5 + + """ + + __doc__ = __doc__.replace("{inherited}", Updater._doc_inherited) + + # Constructor + def __init__(self, trigger, num_swaps, slab_width, target_momentum=math.inf): + # Call the parent constructor with the trigger parameter + super().__init__(trigger) + + # Create a ParameterDict with the given parameters + param_dict = ParameterDict( + num_swaps=int(num_swaps), + slab_width=float(slab_width), + target_momentum=OnlyTypes(float, preprocess=positive_real), + ) + param_dict["target_momentum"] = target_momentum + + # Update the internal parameter dictionary created by the parent class + self._param_dict.update(param_dict) + + # Use the _attach_hook method to create the C++ version of the object + def _attach_hook(self): + sim = self._simulation + self._cpp_obj = _mpcd.ReverseNonequilibriumShearFlow( + sim.state._cpp_sys_def, + self.trigger, + self.num_swaps, + self.slab_width, + self.target_momentum, + ) + + super()._attach_hook() + + @log(category="scalar", requires_run=True) + def summed_exchanged_momentum(self): + r"""float: Total momentum exchanged. + + This quantity logs the total momentum exchanged between all swapped + particle pairs. The value reported is the total exchanged momentum + accumulated till the current timestep. + + """ + return self._cpp_obj.summed_exchanged_momentum + + +__all__ = [ + "ReverseNonequilibriumShearFlow", +] diff --git a/sphinx-doc/hoomd/module-mpcd.rst b/sphinx-doc/hoomd/module-mpcd.rst index c73d197d3b..3a254f6014 100644 --- a/sphinx-doc/hoomd/module-mpcd.rst +++ b/sphinx-doc/hoomd/module-mpcd.rst @@ -17,6 +17,7 @@ mpcd mpcd/module-methods mpcd/module-stream mpcd/module-tune + mpcd/module-update .. rubric:: Classes diff --git a/sphinx-doc/hoomd/mpcd/module-update.rst b/sphinx-doc/hoomd/mpcd/module-update.rst new file mode 100644 index 0000000000..477a38b21e --- /dev/null +++ b/sphinx-doc/hoomd/mpcd/module-update.rst @@ -0,0 +1,13 @@ +update +====== + +.. automodule:: hoomd.mpcd.update + :members: + :exclude-members: ReverseNonequilibriumShearFlow + +.. rubric:: Classes + +.. toctree:: + :maxdepth: 1 + + update/reversenonequilibriumshearflow diff --git a/sphinx-doc/hoomd/mpcd/update/reversenonequilibriumshearflow.rst b/sphinx-doc/hoomd/mpcd/update/reversenonequilibriumshearflow.rst new file mode 100644 index 0000000000..5d13ddc9a0 --- /dev/null +++ b/sphinx-doc/hoomd/mpcd/update/reversenonequilibriumshearflow.rst @@ -0,0 +1,8 @@ +ReverseNonequilibriumShearFlow +============================== + +.. py:currentmodule:: hoomd.mpcd.update + +.. autoclass:: ReverseNonequilibriumShearFlow + :members: + :show-inheritance: diff --git a/sphinx-doc/module-mpcd-update.rst b/sphinx-doc/module-mpcd-update.rst new file mode 100644 index 0000000000..e680e93c9f --- /dev/null +++ b/sphinx-doc/module-mpcd-update.rst @@ -0,0 +1,21 @@ +.. Copyright (c) 2009-2024 The Regents of the University of Michigan. +.. Part of HOOMD-blue, released under the BSD 3-Clause License. + +mpcd.update +----------- + +.. rubric:: Overview + +.. py:currentmodule:: hoomd.mpcd.update + +.. autosummary:: + :nosignatures: + + ReverseNonequilibriumShearFlow + +.. rubric:: Details + +.. automodule:: hoomd.mpcd.update + :synopsis: Updaters. + :members: ReverseNonequilibriumShearFlow + :show-inheritance: From 33a528d74255f7ba658ebb109d8af10d1f0d0ba9 Mon Sep 17 00:00:00 2001 From: Michael Howard Date: Wed, 4 Dec 2024 21:35:48 -0600 Subject: [PATCH 02/10] Implement MPCD reverse nonequilibrium shear flow on GPU --- hoomd/mpcd/CMakeLists.txt | 10 +- hoomd/mpcd/ReverseNonequilibriumShearFlow.cc | 116 +++++----- hoomd/mpcd/ReverseNonequilibriumShearFlow.h | 54 +++-- .../mpcd/ReverseNonequilibriumShearFlowGPU.cc | 207 ++++++++++++++++++ .../mpcd/ReverseNonequilibriumShearFlowGPU.cu | 173 +++++++++++++++ .../ReverseNonequilibriumShearFlowGPU.cuh | 62 ++++++ .../mpcd/ReverseNonequilibriumShearFlowGPU.h | 54 +++++ .../ReverseNonequilibriumShearFlowUtilities.h | 20 +- hoomd/mpcd/module.cc | 2 + hoomd/mpcd/update.py | 90 ++++---- 10 files changed, 655 insertions(+), 133 deletions(-) create mode 100644 hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cc create mode 100644 hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cu create mode 100644 hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cuh create mode 100644 hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.h diff --git a/hoomd/mpcd/CMakeLists.txt b/hoomd/mpcd/CMakeLists.txt index 496a125a95..bf05f7cb5a 100644 --- a/hoomd/mpcd/CMakeLists.txt +++ b/hoomd/mpcd/CMakeLists.txt @@ -76,6 +76,7 @@ if(ENABLE_HIP) CommunicatorGPU.cc ParallelPlateGeometryFillerGPU.cc PlanarPoreGeometryFillerGPU.cc + ReverseNonequilibriumShearFlowGPU.cc SorterGPU.cc SRDCollisionMethodGPU.cc ) @@ -84,6 +85,8 @@ if(ENABLE_HIP) ATCollisionMethodGPU.h BounceBackNVEGPU.cuh BounceBackNVEGPU.h + BounceBackStreamingMethodGPU.cuh + BounceBackStreamingMethodGPU.h BulkStreamingMethodGPU.h CellCommunicator.cuh CellThermoComputeGPU.cuh @@ -92,15 +95,15 @@ if(ENABLE_HIP) CellListGPU.h CommunicatorGPU.cuh CommunicatorGPU.h - BounceBackStreamingMethodGPU.cuh - BounceBackStreamingMethodGPU.h - ParticleData.cuh ParallelPlateGeometryFillerGPU.cuh ParallelPlateGeometryFillerGPU.h + ParticleData.cuh PlanarPoreGeometryFillerGPU.cuh PlanarPoreGeometryFillerGPU.h RejectionVirtualParticleFillerGPU.cuh RejectionVirtualParticleFillerGPU.h + ReverseNonequilibriumShearFlowGPU.cuh + ReverseNonequilibriumShearFlowGPU.h SorterGPU.cuh SorterGPU.h SRDCollisionMethodGPU.cuh @@ -116,6 +119,7 @@ if(ENABLE_HIP) ParallelPlateGeometryFillerGPU.cu PlanarPoreGeometryFillerGPU.cu RejectionVirtualParticleFillerGPU.cu + ReverseNonequilibriumShearFlowGPU.cu SorterGPU.cu SRDCollisionMethodGPU.cu ) diff --git a/hoomd/mpcd/ReverseNonequilibriumShearFlow.cc b/hoomd/mpcd/ReverseNonequilibriumShearFlow.cc index 05130381d0..8c7e523cea 100644 --- a/hoomd/mpcd/ReverseNonequilibriumShearFlow.cc +++ b/hoomd/mpcd/ReverseNonequilibriumShearFlow.cc @@ -26,9 +26,9 @@ mpcd::ReverseNonequilibriumShearFlow::ReverseNonequilibriumShearFlow( Scalar slab_width, Scalar target_momentum) : Updater(sysdef, trigger), m_mpcd_pdata(sysdef->getMPCDParticleData()), m_num_swap(num_swap), - m_slab_width(slab_width), m_summed_momentum_exchange(0), m_num_lo(0), - m_particles_lo(m_exec_conf), m_num_hi(0), m_particles_hi(m_exec_conf), - m_target_momentum(target_momentum), m_update_slabs(true) + m_slab_width(slab_width), m_target_momentum(target_momentum), m_summed_momentum_exchange(0), + m_num_lo(0), m_particles_lo(m_exec_conf), m_num_hi(0), m_particles_hi(m_exec_conf), + m_update_slabs(true) { m_exec_conf->msg->notice(5) << "Constructing ReverseNonequilibriumShearFlow" << std::endl; @@ -102,6 +102,7 @@ void mpcd::ReverseNonequilibriumShearFlow::update(uint64_t timestep) } findSwapParticles(); + sortOutSwapParticles(); stageSwapParticles(); swapParticleMomentum(); } @@ -181,38 +182,50 @@ void mpcd::ReverseNonequilibriumShearFlow::findSwapParticles() } } while (!filled); + } - // sort the arrays - { - ArrayHandle h_particles_lo(m_particles_lo, - access_location::host, - access_mode::readwrite); - ArrayHandle h_particles_hi(m_particles_hi, - access_location::host, - access_mode::readwrite); - ArrayHandle h_tag(m_mpcd_pdata->getTags(), +void mpcd::ReverseNonequilibriumShearFlow::sortOutSwapParticles() + { + ArrayHandle h_particles_lo(m_particles_lo, access_location::host, - access_mode::read); - - if (std::isinf(m_target_momentum)) - { - std::sort(h_particles_lo.data, - h_particles_lo.data + m_num_lo, - detail::MaximumMomentum(h_tag.data)); - std::sort(h_particles_hi.data, - h_particles_hi.data + m_num_hi, - detail::MinimumMomentum(h_tag.data)); - } - else - { - std::sort(h_particles_lo.data, - h_particles_lo.data + m_num_lo, - detail::CompareMomentumToTarget(m_target_momentum, h_tag.data)); - std::sort(h_particles_hi.data, - h_particles_hi.data + m_num_hi, - detail::CompareMomentumToTarget(-m_target_momentum, h_tag.data)); - } + access_mode::readwrite); + ArrayHandle h_particles_hi(m_particles_hi, + access_location::host, + access_mode::readwrite); + ArrayHandle h_tag(m_mpcd_pdata->getTags(), + access_location::host, + access_mode::read); + + const unsigned int num_top_lo = std::min(m_num_swap, m_num_lo); + const unsigned int num_top_hi = std::min(m_num_swap, m_num_hi); + if (std::isinf(m_target_momentum)) + { + std::partial_sort(h_particles_lo.data, + h_particles_lo.data + num_top_lo, + h_particles_lo.data + m_num_lo, + detail::MaximumMomentum(h_tag.data)); + std::partial_sort(h_particles_hi.data, + h_particles_hi.data + num_top_hi, + h_particles_hi.data + m_num_hi, + detail::MinimumMomentum(h_tag.data)); } + else + { + std::partial_sort(h_particles_lo.data, + h_particles_lo.data + num_top_lo, + h_particles_lo.data + m_num_lo, + detail::CompareMomentumToTarget(m_target_momentum, h_tag.data)); + std::partial_sort(h_particles_hi.data, + h_particles_hi.data + num_top_hi, + h_particles_hi.data + m_num_hi, + detail::CompareMomentumToTarget(-m_target_momentum, h_tag.data)); + } + + m_top_particles_lo.resize(num_top_lo); + std::copy(h_particles_lo.data, h_particles_lo.data + num_top_lo, m_top_particles_lo.begin()); + + m_top_particles_hi.resize(num_top_hi); + std::copy(h_particles_hi.data, h_particles_hi.data + num_top_hi, m_top_particles_hi.begin()); } /*! @@ -221,21 +234,22 @@ void mpcd::ReverseNonequilibriumShearFlow::findSwapParticles() void mpcd::ReverseNonequilibriumShearFlow::stageSwapParticles() { // determine number of pairs to swap - unsigned int num_lo_global = m_num_lo; - unsigned int num_hi_global = m_num_hi; + const unsigned int num_top_lo = static_cast(m_top_particles_lo.size()); + const unsigned int num_top_hi = static_cast(m_top_particles_hi.size()); + unsigned int num_top_lo_global = num_top_lo; + unsigned int num_top_hi_global = num_top_hi; #ifdef ENABLE_MPI if (m_sysdef->isDomainDecomposed()) { auto mpi_comm = m_exec_conf->getMPICommunicator(); - MPI_Allreduce(MPI_IN_PLACE, &num_lo_global, 1, MPI_UNSIGNED, MPI_SUM, mpi_comm); - MPI_Allreduce(MPI_IN_PLACE, &num_hi_global, 1, MPI_UNSIGNED, MPI_SUM, mpi_comm); + MPI_Allreduce(MPI_IN_PLACE, &num_top_lo_global, 1, MPI_UNSIGNED, MPI_SUM, mpi_comm); + MPI_Allreduce(MPI_IN_PLACE, &num_top_hi_global, 1, MPI_UNSIGNED, MPI_SUM, mpi_comm); } #endif // ENABLE_MPI - const unsigned int num_pairs = std::min(m_num_swap, std::min(num_lo_global, num_hi_global)); + const unsigned int num_pairs + = std::min(m_num_swap, std::min(num_top_lo_global, num_top_hi_global)); // stage swaps into queue - ArrayHandle h_particles_lo(m_particles_lo, access_location::host, access_mode::read); - ArrayHandle h_particles_hi(m_particles_hi, access_location::host, access_mode::read); ArrayHandle h_particles_staged(m_particles_staged, access_location::host, access_mode::overwrite); @@ -257,13 +271,13 @@ void mpcd::ReverseNonequilibriumShearFlow::stageSwapParticles() if (std::isinf(m_target_momentum)) { lo_op = MPI_MAXLOC; - lo_swap.s = (lo_idx < m_num_lo) ? h_particles_lo.data[lo_idx].y : Scalar(0.0); + lo_swap.s = (lo_idx < num_top_lo) ? m_top_particles_lo[lo_idx].y : Scalar(0.0); } else { lo_op = MPI_MINLOC; - lo_swap.s = (lo_idx < m_num_lo) - ? std::fabs(h_particles_lo.data[lo_idx].y - m_target_momentum) + lo_swap.s = (lo_idx < num_top_lo) + ? std::fabs(m_top_particles_lo[lo_idx].y - m_target_momentum) : std::numeric_limits::max(); } MPI_Allreduce(MPI_IN_PLACE, &lo_swap, 1, MPI_HOOMD_SCALAR_INT, lo_op, mpi_comm); @@ -271,11 +285,11 @@ void mpcd::ReverseNonequilibriumShearFlow::stageSwapParticles() // find particle in hi slab to swap (p < 0) Scalar_Int hi_swap; hi_swap.i = rank; - if (hi_idx < m_num_hi) + if (hi_idx < num_top_hi) { hi_swap.s = (std::isinf(m_target_momentum)) - ? h_particles_hi.data[hi_idx].y - : std::fabs(h_particles_hi.data[hi_idx].y + m_target_momentum); + ? m_top_particles_hi[hi_idx].y + : std::fabs(m_top_particles_hi[hi_idx].y + m_target_momentum); } else { @@ -291,12 +305,12 @@ void mpcd::ReverseNonequilibriumShearFlow::stageSwapParticles() unsigned int dest; if (lo_swap.i == rank) { - particle = h_particles_lo.data[lo_idx++]; + particle = m_top_particles_lo[lo_idx++]; dest = hi_swap.i; } else { - particle = h_particles_hi.data[hi_idx++]; + particle = m_top_particles_hi[hi_idx++]; dest = lo_swap.i; } @@ -314,8 +328,8 @@ void mpcd::ReverseNonequilibriumShearFlow::stageSwapParticles() } else // particles are on the same rank { - Scalar2 lo_particle = h_particles_lo.data[lo_idx++]; - Scalar2 hi_particle = h_particles_hi.data[hi_idx++]; + Scalar2 lo_particle = m_top_particles_lo[lo_idx++]; + Scalar2 hi_particle = m_top_particles_hi[hi_idx++]; std::swap(lo_particle.y, hi_particle.y); h_particles_staged.data[m_num_staged++] = lo_particle; h_particles_staged.data[m_num_staged++] = hi_particle; @@ -325,8 +339,8 @@ void mpcd::ReverseNonequilibriumShearFlow::stageSwapParticles() else #endif // ENABLE_MPI { - Scalar2 lo_particle = h_particles_lo.data[lo_idx++]; - Scalar2 hi_particle = h_particles_hi.data[hi_idx++]; + Scalar2 lo_particle = m_top_particles_lo[lo_idx++]; + Scalar2 hi_particle = m_top_particles_hi[hi_idx++]; std::swap(lo_particle.y, hi_particle.y); h_particles_staged.data[m_num_staged++] = lo_particle; h_particles_staged.data[m_num_staged++] = hi_particle; diff --git a/hoomd/mpcd/ReverseNonequilibriumShearFlow.h b/hoomd/mpcd/ReverseNonequilibriumShearFlow.h index c56743438a..f9adab41ae 100644 --- a/hoomd/mpcd/ReverseNonequilibriumShearFlow.h +++ b/hoomd/mpcd/ReverseNonequilibriumShearFlow.h @@ -2,14 +2,14 @@ // Part of HOOMD-blue, released under the BSD 3-Clause License. /*! - * \file ReverseNonequilibriumShearFlow.h - * \brief Declaration of Reverse nonequilibrium shear flow + * \file mpcd/ReverseNonequilibriumShearFlow.h + * \brief Declaration of reverse nonequilibrium shear flow updater. */ #ifndef MPCD_REVERSE_NONEQUILIBRIUM_SHEAR_FLOW_H_ #define MPCD_REVERSE_NONEQUILIBRIUM_SHEAR_FLOW_H_ -#ifdef NVCC +#ifdef __HIPCC__ #error This header cannot be compiled by nvcc #endif @@ -42,7 +42,7 @@ class PYBIND11_EXPORT ReverseNonequilibriumShearFlow : public Updater //! Apply velocity swaps virtual void update(uint64_t timestep); - //! Get max number of swaps + //! Get the maximum number of swapped pairs Scalar getNumSwap() const { return m_num_swap; @@ -51,7 +51,7 @@ class PYBIND11_EXPORT ReverseNonequilibriumShearFlow : public Updater //! Set the maximum number of swapped pairs void setNumSwap(unsigned int num_swap); - //! Get slab width + //! Get the slab width Scalar getSlabWidth() const { return m_slab_width; @@ -60,7 +60,7 @@ class PYBIND11_EXPORT ReverseNonequilibriumShearFlow : public Updater //! Set the slab width void setSlabWidth(Scalar slab_width); - //! Get target momentum + //! Get the target momentum Scalar getTargetMomentum() const { return m_target_momentum; @@ -77,28 +77,34 @@ class PYBIND11_EXPORT ReverseNonequilibriumShearFlow : public Updater protected: std::shared_ptr m_mpcd_pdata; //!< MPCD particle data - unsigned int m_num_swap; //!< maximum number of swaps - Scalar m_slab_width; //!< width of slabs - Scalar m_summed_momentum_exchange; //!< summed momentum excange between slabs - Scalar2 m_pos_lo; //!< position of bottom slab in box - Scalar2 m_pos_hi; //!< position of top slab in box - unsigned int m_num_lo; //!< number of particles in bottom slab - GPUArray m_particles_lo; //!< List of all particles (indices,momentum) in bottom slab - //!< sorted by momentum closest to +m_target_momentum - unsigned int m_num_hi; //!< number of particles in top slab - GPUArray m_particles_hi; //!< List of all particles (indices,momentum) in top slab - //!< sorted by momentum closest to -m_target_momentum - unsigned int m_num_staged; //!< number of particles staged for swapping - GPUArray m_particles_staged; //!< List of all particles staged for swapping - Scalar m_target_momentum; //!< target momentum for particles in the slabs - - //! Find candidate particles for swapping in the slabs + + unsigned int m_num_swap; //!< Maximum number of swaps + Scalar m_slab_width; //!< Width of slabs + Scalar m_target_momentum; //!< Target momentum + + Scalar m_summed_momentum_exchange; //!< Total momentum exchanged so far + Scalar2 m_pos_lo; //!< y bounds of lower slab + Scalar2 m_pos_hi; //!< y bounds of upper slab + unsigned int m_num_lo; //!< Number of particles in lower slab + GPUArray m_particles_lo; //!< Sorted particle indexes and momenta in lower slab + unsigned int m_num_hi; //!< Number of particles in upper slab + GPUArray m_particles_hi; //!< Sorted particle indexes and momenta in upper slab + + std::vector m_top_particles_lo; //!< Top candidates for swapping in lower slab + std::vector m_top_particles_hi; //!< Top candidates for swapping in upper slab + unsigned int m_num_staged; //!< Number of particles staged for swapping + GPUArray m_particles_staged; //!< Particle indexes and momenta staged for swapping + + //! Find candidate particles for swapping virtual void findSwapParticles(); - //! Stage particle momentum for swapping + //! Sort and copy only the top candidates for swapping + virtual void sortOutSwapParticles(); + + //! Stage particles for swapping void stageSwapParticles(); - //! Swaps momentum between the slabs + //! Swap particle momenta virtual void swapParticleMomentum(); private: diff --git a/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cc b/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cc new file mode 100644 index 0000000000..00ecbc6a88 --- /dev/null +++ b/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cc @@ -0,0 +1,207 @@ +// Copyright (c) 2009-2024 The Regents of the University of Michigan. +// Part of HOOMD-blue, released under the BSD 3-Clause License. + +/*! + * \file mpcd/ReverseNonequilibriumShearFlowGPU.h + * \brief Definition of mpcd::ReverseNonequilibriumShearFlowGPU + */ + +#include "ReverseNonequilibriumShearFlowGPU.cuh" +#include "ReverseNonequilibriumShearFlowGPU.h" +#include "ReverseNonequilibriumShearFlowUtilities.h" + +namespace hoomd + { + +mpcd::ReverseNonequilibriumShearFlowGPU::ReverseNonequilibriumShearFlowGPU( + std::shared_ptr sysdef, + std::shared_ptr trigger, + unsigned int num_swap, + Scalar slab_width, + Scalar target_momentum) + : mpcd::ReverseNonequilibriumShearFlow(sysdef, trigger, num_swap, slab_width, target_momentum), + m_num_lo_hi(2, m_exec_conf), m_momentum_sum(m_exec_conf) + { + m_tuner_filter.reset(new Autotuner<1>({AutotunerBase::makeBlockSizeRange(m_exec_conf)}, + m_exec_conf, + "mpcd_rnes_filter")); + m_tuner_swap.reset(new Autotuner<1>({AutotunerBase::makeBlockSizeRange(m_exec_conf)}, + m_exec_conf, + "mpcd_rnes_swap")); + m_autotuners.insert(m_autotuners.end(), {m_tuner_filter, m_tuner_swap}); + } + +void mpcd::ReverseNonequilibriumShearFlowGPU::findSwapParticles() + { + // fill the layers with (unsorted) particles. this uses do loop for reallocation + bool filled = false; + do + { + const unsigned int num_lo_alloc + = static_cast(m_particles_lo.getNumElements()); + const unsigned int num_hi_alloc + = static_cast(m_particles_hi.getNumElements()); + // filter particles into their slab in y-direction and record momentum in x-direction + { + ArrayHandle d_pos(m_mpcd_pdata->getPositions(), + access_location::device, + access_mode::read); + ArrayHandle d_vel(m_mpcd_pdata->getVelocities(), + access_location::device, + access_mode::read); + + ArrayHandle d_particles_lo(m_particles_lo, + access_location::device, + access_mode::overwrite); + ArrayHandle d_particles_hi(m_particles_hi, + access_location::device, + access_mode::overwrite); + ArrayHandle d_num_lo_hi(m_num_lo_hi, + access_location::device, + access_mode::overwrite); + + m_tuner_filter->begin(); + mpcd::gpu::rnes_filter_particles(d_particles_lo.data, + d_particles_hi.data, + d_num_lo_hi.data, + d_pos.data, + d_vel.data, + m_mpcd_pdata->getMass(), + m_pos_lo, + m_pos_hi, + num_lo_alloc, + num_hi_alloc, + m_mpcd_pdata->getN(), + m_tuner_filter->getParam()[0]); + if (m_exec_conf->isCUDAErrorCheckingEnabled()) + CHECK_CUDA_ERROR(); + m_tuner_filter->end(); + } + + // read number in each slab + { + ArrayHandle h_num_lo_hi(m_num_lo_hi, + access_location::host, + access_mode::read); + m_num_lo = h_num_lo_hi.data[0]; + m_num_hi = h_num_lo_hi.data[1]; + } + + // reallocate if required and go again, this won't happen too much + filled = true; + if (m_num_lo > num_lo_alloc) + { + GPUArray particles_lo(m_num_lo, m_exec_conf); + m_particles_lo.swap(particles_lo); + filled = false; + } + if (m_num_hi > num_hi_alloc) + { + GPUArray particles_hi(m_num_hi, m_exec_conf); + m_particles_hi.swap(particles_hi); + filled = false; + } + } while (!filled); + } + +void mpcd::ReverseNonequilibriumShearFlowGPU::sortOutSwapParticles() + { + ArrayHandle d_particles_lo(m_particles_lo, + access_location::device, + access_mode::readwrite); + ArrayHandle d_particles_hi(m_particles_hi, + access_location::device, + access_mode::readwrite); + ArrayHandle d_tag(m_mpcd_pdata->getTags(), + access_location::device, + access_mode::read); + + if (std::isinf(m_target_momentum)) + { + mpcd::gpu::rnes_sort(d_particles_lo.data, m_num_lo, detail::MaximumMomentum(d_tag.data)); + mpcd::gpu::rnes_sort(d_particles_hi.data, m_num_hi, detail::MinimumMomentum(d_tag.data)); + } + else + { + mpcd::gpu::rnes_sort(d_particles_lo.data, + m_num_lo, + detail::CompareMomentumToTarget(m_target_momentum, d_tag.data)); + mpcd::gpu::rnes_sort(d_particles_hi.data, + m_num_hi, + detail::CompareMomentumToTarget(-m_target_momentum, d_tag.data)); + } + + // copy only the top particles to the host memory + const unsigned int num_top_lo = std::min(m_num_swap, m_num_lo); + m_top_particles_lo.resize(num_top_lo); + mpcd::gpu::rnes_copy_top_particles(m_top_particles_lo.data(), d_particles_lo.data, num_top_lo); + if (m_exec_conf->isCUDAErrorCheckingEnabled()) + CHECK_CUDA_ERROR(); + + const unsigned int num_top_hi = std::min(m_num_swap, m_num_hi); + m_top_particles_hi.resize(num_top_hi); + mpcd::gpu::rnes_copy_top_particles(m_top_particles_hi.data(), d_particles_hi.data, num_top_hi); + if (m_exec_conf->isCUDAErrorCheckingEnabled()) + CHECK_CUDA_ERROR(); + } + +void mpcd::ReverseNonequilibriumShearFlowGPU::swapParticleMomentum() + { + ArrayHandle d_vel(m_mpcd_pdata->getVelocities(), + access_location::device, + access_mode::readwrite); + ArrayHandle d_particles_staged(m_particles_staged, + access_location::device, + access_mode::read); + + // perform swap + m_momentum_sum.resetFlags(0); + m_tuner_swap->begin(); + mpcd::gpu::rnes_swap_particles(d_vel.data, + m_momentum_sum.getDeviceFlags(), + d_particles_staged.data, + m_mpcd_pdata->getMass(), + m_num_staged, + m_tuner_swap->getParam()[0]); + if (m_exec_conf->isCUDAErrorCheckingEnabled()) + CHECK_CUDA_ERROR(); + m_tuner_swap->end(); + + // sum momentum exchange from this swap + Scalar momentum_sum = m_momentum_sum.readFlags(); +#ifdef ENABLE_MPI + if (m_sysdef->isDomainDecomposed()) + { + auto comm = m_exec_conf->getMPICommunicator(); + MPI_Allreduce(MPI_IN_PLACE, &momentum_sum, 1, MPI_HOOMD_SCALAR, MPI_SUM, comm); + } +#endif + + // absolute value gives extra factor of 2, divide out when accumulating + m_summed_momentum_exchange += Scalar(0.5) * momentum_sum; + } + +namespace mpcd + { +namespace detail + { +/*! + * \param m Python module to export to + */ +void export_ReverseNonequilibriumShearFlowGPU(pybind11::module& m) + { + pybind11::class_>( + m, + "ReverseNonequilibriumShearFlowGPU") + .def(pybind11::init, + std::shared_ptr, + unsigned int, + Scalar, + Scalar>()); + } + } // namespace detail + } // namespace mpcd + + } // end namespace hoomd diff --git a/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cu b/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cu new file mode 100644 index 0000000000..1d80310494 --- /dev/null +++ b/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cu @@ -0,0 +1,173 @@ +// Copyright (c) 2009-2024 The Regents of the University of Michigan. +// Part of HOOMD-blue, released under the BSD 3-Clause License. + +/*! + * \file mpcd/ReverseNonequilibriumShearFlowGPU.cu + * \brief Defines GPU functions and kernels used by mpcd::ReverseNonequilibriumShearFlowGPU + */ + +#include "ReverseNonequilibriumShearFlowGPU.cuh" +#include "ReverseNonequilibriumShearFlowUtilities.h" + +namespace hoomd + { +namespace mpcd + { +namespace gpu + { +namespace kernel + { + +__global__ void rnes_filter_particles(Scalar2* d_particles_lo, + Scalar2* d_particles_hi, + unsigned int* d_num_lo_hi, + const Scalar4* d_pos, + const Scalar4* d_vel, + const Scalar mass, + const Scalar2 pos_lo, + const Scalar2 pos_hi, + const unsigned int num_lo_alloc, + const unsigned int num_hi_alloc, + const unsigned int N) + { + // one thread per particle + unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x; + if (idx >= N) + return; + + const Scalar4 vel = d_vel[idx]; + const Scalar momentum = vel.x * mass; + const Scalar y = d_pos[idx].y; + + if (pos_lo.x <= y && y < pos_lo.y + && momentum > Scalar(0.0)) // lower slab, search for positive momentum + { + const unsigned int num_lo = atomicInc(d_num_lo_hi, 0xffffffff); + if (num_lo < num_lo_alloc) + { + d_particles_lo[num_lo] = make_scalar2(__int_as_scalar(idx), momentum); + } + } + else if (pos_hi.x <= y && y < pos_hi.y + && momentum < Scalar(0.0)) // higher slab, search for negative momentum + { + const unsigned int num_hi = atomicInc(d_num_lo_hi + 1, 0xffffffff); + if (num_hi < num_hi_alloc) + { + d_particles_hi[num_hi] = make_scalar2(__int_as_scalar(idx), momentum); + } + } + } + +__global__ void rnes_swap_particles(Scalar4* d_vel, + Scalar* d_momentum_sum, + const Scalar2* d_particles_staged, + const Scalar mass, + const unsigned int num_staged) + { + // one thread per particle + unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x; + if (idx >= num_staged) + return; + + const Scalar2 pidx_mom = d_particles_staged[idx]; + const unsigned int pidx = __scalar_as_int(pidx_mom.x); + const Scalar new_momentum = pidx_mom.y; + + const Scalar current_momentum = d_vel[pidx].x * mass; + + d_vel[pidx].x = new_momentum / mass; + atomicAdd(d_momentum_sum, std::fabs(new_momentum - current_momentum)); + } + + } // end namespace kernel + +cudaError_t rnes_filter_particles(Scalar2* d_particles_lo, + Scalar2* d_particles_hi, + unsigned int* d_num_lo_hi, + const Scalar4* d_pos, + const Scalar4* d_vel, + const Scalar mass, + const Scalar2& pos_lo, + const Scalar2& pos_hi, + const unsigned int num_lo_alloc, + const unsigned int num_hi_alloc, + const unsigned int N, + const unsigned int block_size) + { + // zero counts in each bin + cudaError_t error = cudaMemset(d_num_lo_hi, 0, 2 * sizeof(unsigned int)); + if (error != cudaSuccess) + return error; + + unsigned int max_block_size; + cudaFuncAttributes attr; + cudaFuncGetAttributes(&attr, (const void*)mpcd::gpu::kernel::rnes_filter_particles); + max_block_size = attr.maxThreadsPerBlock; + + unsigned int run_block_size = min(block_size, max_block_size); + + dim3 grid(N / run_block_size + 1); + mpcd::gpu::kernel::rnes_filter_particles<<>>(d_particles_lo, + d_particles_hi, + d_num_lo_hi, + d_pos, + d_vel, + mass, + pos_lo, + pos_hi, + num_lo_alloc, + num_hi_alloc, + N); + + return cudaSuccess; + } + +template void rnes_sort(Scalar2*, + const unsigned int, + const mpcd::detail::MinimumMomentum&); +template void rnes_sort(Scalar2*, + const unsigned int, + const mpcd::detail::MaximumMomentum&); +template void +rnes_sort(Scalar2*, + const unsigned int, + const mpcd::detail::CompareMomentumToTarget&); + +cudaError_t rnes_copy_top_particles(Scalar2* h_top_particles, + const Scalar2* d_particles, + const unsigned int num_top) + { + return cudaMemcpy(h_top_particles, + d_particles, + num_top * sizeof(Scalar2), + cudaMemcpyDeviceToHost); + } + +cudaError_t rnes_swap_particles(Scalar4* d_vel, + Scalar* d_momentum_sum, + const Scalar2* d_particles_staged, + const Scalar mass, + const unsigned int num_staged, + const unsigned int block_size) + { + unsigned int max_block_size; + cudaFuncAttributes attr; + cudaFuncGetAttributes(&attr, (const void*)mpcd::gpu::kernel::rnes_swap_particles); + max_block_size = attr.maxThreadsPerBlock; + + unsigned int run_block_size = min(block_size, max_block_size); + + dim3 grid(num_staged / run_block_size + 1); + mpcd::gpu::kernel::rnes_swap_particles<<>>(d_vel, + d_momentum_sum, + d_particles_staged, + mass, + num_staged); + + return cudaSuccess; + } + + } // end namespace gpu + } // end namespace mpcd + } // end namespace hoomd diff --git a/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cuh b/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cuh new file mode 100644 index 0000000000..627b123852 --- /dev/null +++ b/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cuh @@ -0,0 +1,62 @@ +// Copyright (c) 2009-2024 The Regents of the University of Michigan. +// Part of HOOMD-blue, released under the BSD 3-Clause License. + +#ifndef MPCD_REVERSE_NONEQUILIBRIUM_SHEAR_FLOW_GPU_CUH_ +#define MPCD_REVERSE_NONEQUILIBRIUM_SHEAR_FLOW_GPU_CUH_ + +/*! + * \file mpcd/ReverseNonequilibriumShearFlowGPU.cuh + * \brief Declaration of CUDA kernels for mpcd::ReverseNonequilibriumShearFlowGPU + */ + +#include +#include +#include + +#include "hoomd/HOOMDMath.h" + +namespace hoomd + { +namespace mpcd + { +namespace gpu + { + +//! Filter particles into slabs +cudaError_t rnes_filter_particles(Scalar2* d_particles_lo, + Scalar2* d_particles_hi, + unsigned int* d_num_lo_hi, + const Scalar4* d_pos, + const Scalar4* d_vel, + const Scalar mass, + const Scalar2& pos_lo, + const Scalar2& pos_hi, + const unsigned int num_lo_alloc, + const unsigned int num_hi_alloc, + const unsigned int N, + const unsigned int block_size); + +template void rnes_sort(Scalar2* d_particles, const unsigned int N, const T& comp); + +cudaError_t rnes_copy_top_particles(Scalar2* h_top_particles, + const Scalar2* d_particles, + const unsigned int num_top); + +cudaError_t rnes_swap_particles(Scalar4* d_vel, + Scalar* d_momentum_sum, + const Scalar2* d_particles_staged, + const Scalar mass, + const unsigned int num_staged, + const unsigned int block_size); + +#ifdef __HIPCC__ +template void rnes_sort(Scalar2* d_particles, const unsigned int N, const T& comp) + { + thrust::sort(thrust::device, d_particles, d_particles + N, comp); + } +#endif + + } // end namespace gpu + } // end namespace mpcd + } // end namespace hoomd +#endif // MPCD_REVERSE_NONEQUILIBRIUM_SHEAR_FLOW_GPU_CUH_ diff --git a/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.h b/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.h new file mode 100644 index 0000000000..8b7c930fc3 --- /dev/null +++ b/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.h @@ -0,0 +1,54 @@ +// Copyright (c) 2009-2024 The Regents of the University of Michigan. +// Part of HOOMD-blue, released under the BSD 3-Clause License. + +/*! + * \file mpcd/ReverseNonequilibriumShearFlowGPU.h + * \brief Declaration of reverse nonequilibrium shear flow updater on the GPU. + */ + +#ifndef MPCD_REVERSE_NONEQUILIBRIUM_SHEAR_FLOW_GPU_H_ +#define MPCD_REVERSE_NONEQUILIBRIUM_SHEAR_FLOW_GPU_H_ + +#ifdef __HIPCC__ +#error This header cannot be compiled by nvcc +#endif + +#include "ReverseNonequilibriumShearFlow.h" +#include "hoomd/Autotuner.h" +#include "hoomd/GPUFlags.h" +#include + +namespace hoomd + { +namespace mpcd + { +//! Reverse nonequilibrium shear flow updater on GPU +class PYBIND11_EXPORT ReverseNonequilibriumShearFlowGPU : public ReverseNonequilibriumShearFlow + { + public: + //! Constructor + ReverseNonequilibriumShearFlowGPU(std::shared_ptr sysdef, + std::shared_ptr trigger, + unsigned int num_swap, + Scalar slab_width, + Scalar target_momentum); + + protected: + //! Find candidate particles for swapping + void findSwapParticles() override; + + //! Sort and copy only the top candidates for swapping + void sortOutSwapParticles() override; + + //! Swap particle momenta + void swapParticleMomentum() override; + + private: + GPUArray m_num_lo_hi; + GPUFlags m_momentum_sum; + std::shared_ptr> m_tuner_filter; //!< Autotuner for filtering particles + std::shared_ptr> m_tuner_swap; //!< Autotuner for swapping particles + }; + } // end namespace mpcd + } // end namespace hoomd +#endif // MPCD_REVERSE_NONEQUILIBRIUM_SHEAR_FLOW_GPU_H_ diff --git a/hoomd/mpcd/ReverseNonequilibriumShearFlowUtilities.h b/hoomd/mpcd/ReverseNonequilibriumShearFlowUtilities.h index f88a8184cb..9cfbe04c9b 100644 --- a/hoomd/mpcd/ReverseNonequilibriumShearFlowUtilities.h +++ b/hoomd/mpcd/ReverseNonequilibriumShearFlowUtilities.h @@ -11,6 +11,12 @@ #include "hoomd/HOOMDMath.h" +#ifdef __HIPCC__ +#define HOSTDEVICE __host__ __device__ inline +#else +#define HOSTDEVICE inline __attribute__((always_inline)) +#endif // __HIPCC__ + namespace hoomd { namespace mpcd @@ -21,12 +27,12 @@ namespace detail class CompareMomentumToTarget { public: - CompareMomentumToTarget(Scalar target_momentum, const unsigned int* tags_) + HOSTDEVICE CompareMomentumToTarget(Scalar target_momentum, const unsigned int* tags_) : p(target_momentum), tags(tags_) { } - bool operator()(const Scalar2& in0, const Scalar2& in1) const + HOSTDEVICE bool operator()(const Scalar2& in0, const Scalar2& in1) const { const Scalar dp0 = std::fabs(in0.y - p); const Scalar dp1 = std::fabs(in1.y - p); @@ -55,9 +61,9 @@ class CompareMomentumToTarget class MaximumMomentum { public: - MaximumMomentum(const unsigned int* tags_) : tags(tags_) { } + HOSTDEVICE MaximumMomentum(const unsigned int* tags_) : tags(tags_) { } - bool operator()(const Scalar2& in0, const Scalar2& in1) const + HOSTDEVICE bool operator()(const Scalar2& in0, const Scalar2& in1) const { const Scalar p0 = in0.y; const Scalar p1 = in1.y; @@ -85,9 +91,9 @@ class MaximumMomentum class MinimumMomentum { public: - MinimumMomentum(const unsigned int* tags_) : tags(tags_) { } + HOSTDEVICE MinimumMomentum(const unsigned int* tags_) : tags(tags_) { } - bool operator()(const Scalar2& in0, const Scalar2& in1) const + HOSTDEVICE bool operator()(const Scalar2& in0, const Scalar2& in1) const { const Scalar p0 = in0.y; const Scalar p1 = in1.y; @@ -116,4 +122,6 @@ class MinimumMomentum } // end namespace mpcd } // end namespace hoomd +#undef HOSTDEVICE + #endif // MPCD_REVERSE_NONEQUILIBRIUM_SHEAR_FLOW_UTILITIES_H_ diff --git a/hoomd/mpcd/module.cc b/hoomd/mpcd/module.cc index d97d8bb400..d765cca3ed 100644 --- a/hoomd/mpcd/module.cc +++ b/hoomd/mpcd/module.cc @@ -62,6 +62,7 @@ void export_CosineChannelGeometryFillerGPU(pybind11::module&); void export_CosineExpansionContractionGeometryFillerGPU(pybind11::module&); void export_ParallelPlateGeometryFillerGPU(pybind11::module&); void export_PlanarPoreGeometryFillerGPU(pybind11::module&); +void export_ReverseNonequilibriumShearFlowGPU(pybind11::module&); void export_SorterGPU(pybind11::module&); void export_SphereGeometryFillerGPU(pybind11::module&); void export_SRDCollisionMethodGPU(pybind11::module&); @@ -237,6 +238,7 @@ PYBIND11_MODULE(_mpcd, m) export_CosineExpansionContractionGeometryFillerGPU(m); export_ParallelPlateGeometryFillerGPU(m); export_PlanarPoreGeometryFillerGPU(m); + export_ReverseNonequilibriumShearFlowGPU(m); export_SorterGPU(m); export_SphereGeometryFillerGPU(m); export_SRDCollisionMethodGPU(m); diff --git a/hoomd/mpcd/update.py b/hoomd/mpcd/update.py index 296866ee40..960d02537a 100644 --- a/hoomd/mpcd/update.py +++ b/hoomd/mpcd/update.py @@ -10,7 +10,8 @@ """ -from . import _mpcd +import hoomd +from hoomd.mpcd import _mpcd from hoomd.operation import Updater from hoomd.data.parameterdicts import ParameterDict from hoomd.logging import log @@ -20,49 +21,41 @@ class ReverseNonequilibriumShearFlow(Updater): - r"""MPCD reverse nonequilibrium shear flow. + r"""Reverse nonequilibrium shear flow. Args: - trigger (hoomd.trigger.trigger_like): Trigger to swap momentum. + trigger (hoomd.trigger.trigger_like): Select the time steps on which to + to swap momentum. - num_swaps (int): Maximum number of times to swap momentum per update. + num_swaps (int): Maximum number of pairs to swap per update. slab_width (float): Width of momentum-exchange slabs. - target_momentum (float): Target momentum for swapped particles. This - argument has a default value of infinity but can be - redefined with positive real numbers only. - - This updater generates a bidirectional shear flow in *x* by imposing a - momentum flux on the system in *y*. - There are two exchange slabs separated by a distance of :math:`L_y/2`. The - edges of these slabs are located at (:math:`-L_y/2`, :math:`-L_y/2` + - `slab_width`) and (:math:`0.0`, `slab_width`) along the *y*-direction - (gradient direction) of the simulation box. On each `trigger`, particles - are sorted into the exchange slabs based on their positions in the box. - Particles whose *x*-component momenta are near `target_momentum` - in the lower slab, and those near -`target_momentum` in the upper - slab, are selected for a pairwise momentum swap. Up to `num_swaps` swaps - are executed per update. - - The amount of momentum transferred from the lower slab to the upper slab is - known. Therefore `summed_exchanged_momentum`, which returns the accumulated - momentum exchanged till the current timestep, is used to calculate the - momentum flux, :math:`j_{xy}`. The shear rate, :math:`\dot{\gamma}`, is - also extracted as the gradient of the linear velocity profile developed - from the flow. The viscosity can then be computed from these two quantities - as: - - .. math:: - - \eta (\dot{\gamma}) = \frac{j_{xy}}{\dot{\gamma}}. + target_momentum (float): Target magnitude of momentum for swapped + particles (must be positive). + + `ReverseNonequilibriumShearFlow` generates a bidirectional shear flow in + *x* by imposing a momentum flux on MPCD particles in *y*. Particles are + selected from two momentum-exchange slabs with normal to *y*, width *w*, + and separated by :math:`L_y/2`. The lower slab accordingly has particles + with :math:`-L_y/2 \le y < L_y/2 + w`, while the upper slab has particles + with :math:`0 \le y < w`. + + MPCD particles are sorted into these slabs, and the particles whose *x* + momenta are closest to the `target_momentum` :math:`p_0` in the lower slab + and :math:`-p_0` in the upper slab are selected for a momentum swap. + Up to `num_swaps` swaps are executed each time. + + The amount of momentum transferred from the lower slab to the upper slab + is accumulated into `summed_exchangeD_momentum`. This quantity can be used + to calculate the momentum flux and, in conjunction with the shear velocity + field that is generated, determine the shear viscosity. .. rubric:: Examples: - In the original implementation by - `Müller-Plathe `_, - only the fastest particle and the slowest particle are swapped. This is - achieved by setting `num_swaps` to 1 and keeping `target_momentum` at its + To implement the method as originally proposed by `Müller-Plathe`_, + only the fastest particle and the slowest particle in the momentum-exchange + slabs are swapped. Set `num_swaps` to 1 and `target_momentum` at its default value of infinity. .. code-block:: python @@ -72,8 +65,7 @@ class ReverseNonequilibriumShearFlow(Updater): ) simulation.operations.updaters.append(flow) - An alternative approach proposed by - `Tenney and Maginn `_ swaps particles + An alternative approach proposed by `Tenney and Maginn`_ swaps particles that are instead closest to the `target_momentum`, typically requiring more swaps per update. @@ -123,31 +115,32 @@ class ReverseNonequilibriumShearFlow(Updater): flow.target_momentum = 5 + .. _Müller-Plathe: https://doi.org/10.1103/PhysRevE.59.4894 + .. _Tenney and Maginn: https://doi.org/10.1063/1.3276454 + """ __doc__ = __doc__.replace("{inherited}", Updater._doc_inherited) - # Constructor def __init__(self, trigger, num_swaps, slab_width, target_momentum=math.inf): - # Call the parent constructor with the trigger parameter super().__init__(trigger) - # Create a ParameterDict with the given parameters param_dict = ParameterDict( num_swaps=int(num_swaps), slab_width=float(slab_width), target_momentum=OnlyTypes(float, preprocess=positive_real), ) param_dict["target_momentum"] = target_momentum - - # Update the internal parameter dictionary created by the parent class self._param_dict.update(param_dict) - # Use the _attach_hook method to create the C++ version of the object def _attach_hook(self): - sim = self._simulation - self._cpp_obj = _mpcd.ReverseNonequilibriumShearFlow( - sim.state._cpp_sys_def, + if isinstance(self._simulation.device, hoomd.device.GPU): + class_ = _mpcd.ReverseNonequilibriumShearFlowGPU + else: + class_ = _mpcd.ReverseNonequilibriumShearFlow + + self._cpp_obj = class_( + self._simulation.state._cpp_sys_def, self.trigger, self.num_swaps, self.slab_width, @@ -160,9 +153,8 @@ def _attach_hook(self): def summed_exchanged_momentum(self): r"""float: Total momentum exchanged. - This quantity logs the total momentum exchanged between all swapped - particle pairs. The value reported is the total exchanged momentum - accumulated till the current timestep. + This quantity is the total momentum exchanged from the lower slab + to the upper slab. """ return self._cpp_obj.summed_exchanged_momentum From 869043ca144e291d6fb5730eb79f627cf1088435 Mon Sep 17 00:00:00 2001 From: Michael Howard Date: Thu, 12 Dec 2024 20:20:36 -0600 Subject: [PATCH 03/10] Only include thrust under nvcc --- hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cuh | 2 ++ 1 file changed, 2 insertions(+) diff --git a/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cuh b/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cuh index 627b123852..6924634ab3 100644 --- a/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cuh +++ b/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cuh @@ -10,8 +10,10 @@ */ #include +#ifdef __HIPCC__ #include #include +#endif // __HIPCC__ #include "hoomd/HOOMDMath.h" From 4714b00f63e9976f3a41fff56dd0207722308c7f Mon Sep 17 00:00:00 2001 From: Michael Howard Date: Thu, 2 Jan 2025 16:20:05 -0600 Subject: [PATCH 04/10] Move thrust code to cu file --- .../mpcd/ReverseNonequilibriumShearFlowGPU.cu | 40 ++++++++++++++----- .../ReverseNonequilibriumShearFlowGPU.cuh | 17 ++++---- 2 files changed, 37 insertions(+), 20 deletions(-) diff --git a/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cu b/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cu index 1d80310494..add25b12f8 100644 --- a/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cu +++ b/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cu @@ -7,7 +7,9 @@ */ #include "ReverseNonequilibriumShearFlowGPU.cuh" -#include "ReverseNonequilibriumShearFlowUtilities.h" + +#include +#include namespace hoomd { @@ -123,16 +125,32 @@ cudaError_t rnes_filter_particles(Scalar2* d_particles_lo, return cudaSuccess; } -template void rnes_sort(Scalar2*, - const unsigned int, - const mpcd::detail::MinimumMomentum&); -template void rnes_sort(Scalar2*, - const unsigned int, - const mpcd::detail::MaximumMomentum&); -template void -rnes_sort(Scalar2*, - const unsigned int, - const mpcd::detail::CompareMomentumToTarget&); +/* + * thrust doesn't like to be included in host code, so explicitly list the ones + * we need as wrappers to a template. + */ +template void rnes_thrust_sort(Scalar2* d_particles, const unsigned int N, const T& comp) + { + thrust::sort(thrust::device, d_particles, d_particles + N, comp); + } +void rnes_sort(Scalar2* d_particles, + const unsigned int N, + const mpcd::detail::MinimumMomentum& comp) + { + rnes_thrust_sort(d_particles, N, comp); + } +void rnes_sort(Scalar2* d_particles, + const unsigned int N, + const mpcd::detail::MaximumMomentum& comp) + { + rnes_thrust_sort(d_particles, N, comp); + } +void rnes_sort(Scalar2* d_particles, + const unsigned int N, + const mpcd::detail::CompareMomentumToTarget& comp) + { + rnes_thrust_sort(d_particles, N, comp); + } cudaError_t rnes_copy_top_particles(Scalar2* h_top_particles, const Scalar2* d_particles, diff --git a/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cuh b/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cuh index 6924634ab3..681a51e514 100644 --- a/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cuh +++ b/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cuh @@ -10,13 +10,11 @@ */ #include -#ifdef __HIPCC__ -#include -#include -#endif // __HIPCC__ #include "hoomd/HOOMDMath.h" +#include "ReverseNonequilibriumShearFlowUtilities.h" + namespace hoomd { namespace mpcd @@ -38,7 +36,11 @@ cudaError_t rnes_filter_particles(Scalar2* d_particles_lo, const unsigned int N, const unsigned int block_size); -template void rnes_sort(Scalar2* d_particles, const unsigned int N, const T& comp); +void rnes_sort(Scalar2* d_particles, const unsigned int N, const detail::MinimumMomentum& comp); +void rnes_sort(Scalar2* d_particles, const unsigned int N, const detail::MaximumMomentum& comp); +void rnes_sort(Scalar2* d_particles, + const unsigned int N, + const detail::CompareMomentumToTarget& comp); cudaError_t rnes_copy_top_particles(Scalar2* h_top_particles, const Scalar2* d_particles, @@ -52,10 +54,7 @@ cudaError_t rnes_swap_particles(Scalar4* d_vel, const unsigned int block_size); #ifdef __HIPCC__ -template void rnes_sort(Scalar2* d_particles, const unsigned int N, const T& comp) - { - thrust::sort(thrust::device, d_particles, d_particles + N, comp); - } + #endif } // end namespace gpu From 912ac3e23d4405385118c8bd0903b49b216af2dc Mon Sep 17 00:00:00 2001 From: "W. Kwabena Darko" <84714936+wkdarko@users.noreply.github.com> Date: Mon, 13 Jan 2025 14:27:48 -0600 Subject: [PATCH 05/10] corrected typo in update.py --- hoomd/mpcd/update.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hoomd/mpcd/update.py b/hoomd/mpcd/update.py index 960d02537a..7485728c76 100644 --- a/hoomd/mpcd/update.py +++ b/hoomd/mpcd/update.py @@ -47,7 +47,7 @@ class ReverseNonequilibriumShearFlow(Updater): Up to `num_swaps` swaps are executed each time. The amount of momentum transferred from the lower slab to the upper slab - is accumulated into `summed_exchangeD_momentum`. This quantity can be used + is accumulated into `summed_exchanged_momentum`. This quantity can be used to calculate the momentum flux and, in conjunction with the shear velocity field that is generated, determine the shear viscosity. From 206d5d20513cc7586cc2c692ac7519a7aea7389a Mon Sep 17 00:00:00 2001 From: "W. Kwabena Darko" <84714936+wkdarko@users.noreply.github.com> Date: Tue, 14 Jan 2025 08:29:34 -0600 Subject: [PATCH 06/10] Update credits.rst --- sphinx-doc/credits.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/sphinx-doc/credits.rst b/sphinx-doc/credits.rst index b841388e32..cd9f60e9c1 100644 --- a/sphinx-doc/credits.rst +++ b/sphinx-doc/credits.rst @@ -72,6 +72,7 @@ The following people have contributed to HOOMD-blue: * Khalid Ahmed, University of Michigan * Kody Takada, University of Michigan * Kristi Pepa, University of Michigan +* Kwabena Darko, University of Houston * Kwanghwi Je, University of Michigan * Kieran Nehil-Puleo, Vanderbilt University * Lin Yang, Iowa State University From 089f8a03652e6270565ed244d36751028925ee14 Mon Sep 17 00:00:00 2001 From: "W. Kwabena Darko" <84714936+wkdarko@users.noreply.github.com> Date: Tue, 14 Jan 2025 09:14:12 -0600 Subject: [PATCH 07/10] Update CHANGELOG.rst --- CHANGELOG.rst | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 66851357f7..79dc0f658b 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -4,6 +4,15 @@ Change Log 5.x --- +5.0.1 (not yet released) +^^^^^^^^^^^^^^^^^^^^^^^^ + +*Added* + +* ``mpcd.update.ReverseNonequilibriumShearFlow`` + (`#1983 `__). + + 5.0.0 (2024-12-02) ^^^^^^^^^^^^^^^^^^ From 2c6329205e3526fcc20378efeee677f9d81aed81 Mon Sep 17 00:00:00 2001 From: "W. Kwabena Darko" <84714936+wkdarko@users.noreply.github.com> Date: Tue, 14 Jan 2025 11:14:36 -0600 Subject: [PATCH 08/10] Apply suggestions from code review Co-authored-by: Michael Howard --- CHANGELOG.rst | 2 +- hoomd/mpcd/ReverseNonequilibriumShearFlow.cc | 6 ------ 2 files changed, 1 insertion(+), 7 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 79dc0f658b..8b14329128 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -4,7 +4,7 @@ Change Log 5.x --- -5.0.1 (not yet released) +5.1.0 (not yet released) ^^^^^^^^^^^^^^^^^^^^^^^^ *Added* diff --git a/hoomd/mpcd/ReverseNonequilibriumShearFlow.cc b/hoomd/mpcd/ReverseNonequilibriumShearFlow.cc index 8c7e523cea..908c76d04d 100644 --- a/hoomd/mpcd/ReverseNonequilibriumShearFlow.cc +++ b/hoomd/mpcd/ReverseNonequilibriumShearFlow.cc @@ -13,12 +13,6 @@ namespace hoomd { -/*! - * \param sysdata MPCD system data this updater will act on - * \param num_swap Max number of swaps - * \param slab_width Slab width - * \param target_momentum Target momentum in x-direction for the two slabs - */ mpcd::ReverseNonequilibriumShearFlow::ReverseNonequilibriumShearFlow( std::shared_ptr sysdef, std::shared_ptr trigger, From 210a0c09c520bf1fd3db7357401c34d30eff30d5 Mon Sep 17 00:00:00 2001 From: "W. Kwabena Darko" Date: Tue, 14 Jan 2025 13:40:16 -0600 Subject: [PATCH 09/10] slab check generalized for triclinic boxes --- hoomd/mpcd/ReverseNonequilibriumShearFlow.cc | 2 +- hoomd/mpcd/update.py | 8 -------- sphinx-doc/module-mpcd-update.rst | 21 -------------------- 3 files changed, 1 insertion(+), 30 deletions(-) delete mode 100644 sphinx-doc/module-mpcd-update.rst diff --git a/hoomd/mpcd/ReverseNonequilibriumShearFlow.cc b/hoomd/mpcd/ReverseNonequilibriumShearFlow.cc index 908c76d04d..d05a52021f 100644 --- a/hoomd/mpcd/ReverseNonequilibriumShearFlow.cc +++ b/hoomd/mpcd/ReverseNonequilibriumShearFlow.cc @@ -73,7 +73,7 @@ void mpcd::ReverseNonequilibriumShearFlow::setSlabWidth(Scalar slab_width) void mpcd::ReverseNonequilibriumShearFlow::setSlabs() { const BoxDim& global_box = m_pdata->getGlobalBox(); - if (m_slab_width > Scalar(0.5) * global_box.getL().y) + if (m_slab_width > Scalar(0.5) * global_box.getNearestPlaneDistance().y) { throw std::runtime_error("Slab width cannot be larger than Ly/2"); } diff --git a/hoomd/mpcd/update.py b/hoomd/mpcd/update.py index 7485728c76..43e18c36ee 100644 --- a/hoomd/mpcd/update.py +++ b/hoomd/mpcd/update.py @@ -83,14 +83,6 @@ class ReverseNonequilibriumShearFlow(Updater): **Members defined in** `ReverseNonequilibriumShearFlow`: Attributes: - trigger (hoomd.trigger.trigger_like): Trigger to swap momentum. - - .. rubric:: Example: - - .. code-block:: python - - flow.trigger = 1 - num_swaps (int): Maximum number of times to swap momentum per update. .. rubric:: Example: diff --git a/sphinx-doc/module-mpcd-update.rst b/sphinx-doc/module-mpcd-update.rst deleted file mode 100644 index e680e93c9f..0000000000 --- a/sphinx-doc/module-mpcd-update.rst +++ /dev/null @@ -1,21 +0,0 @@ -.. Copyright (c) 2009-2024 The Regents of the University of Michigan. -.. Part of HOOMD-blue, released under the BSD 3-Clause License. - -mpcd.update ------------ - -.. rubric:: Overview - -.. py:currentmodule:: hoomd.mpcd.update - -.. autosummary:: - :nosignatures: - - ReverseNonequilibriumShearFlow - -.. rubric:: Details - -.. automodule:: hoomd.mpcd.update - :synopsis: Updaters. - :members: ReverseNonequilibriumShearFlow - :show-inheritance: From f7b6f97fb95ce657104c967cbcab6b3f998bddab Mon Sep 17 00:00:00 2001 From: "W. Kwabena Darko" Date: Mon, 20 Jan 2025 13:23:38 -0600 Subject: [PATCH 10/10] update license headers --- hoomd/mpcd/ReverseNonequilibriumShearFlow.cc | 2 +- hoomd/mpcd/ReverseNonequilibriumShearFlow.h | 2 +- hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cc | 2 +- hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cu | 2 +- hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cuh | 2 +- hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.h | 2 +- hoomd/mpcd/ReverseNonequilibriumShearFlowUtilities.h | 2 +- hoomd/mpcd/pytest/test_update.py | 2 +- hoomd/mpcd/update.py | 2 +- 9 files changed, 9 insertions(+), 9 deletions(-) diff --git a/hoomd/mpcd/ReverseNonequilibriumShearFlow.cc b/hoomd/mpcd/ReverseNonequilibriumShearFlow.cc index d05a52021f..860882eeec 100644 --- a/hoomd/mpcd/ReverseNonequilibriumShearFlow.cc +++ b/hoomd/mpcd/ReverseNonequilibriumShearFlow.cc @@ -1,4 +1,4 @@ -// Copyright (c) 2009-2024 The Regents of the University of Michigan. +// Copyright (c) 2009-2025 The Regents of the University of Michigan. // Part of HOOMD-blue, released under the BSD 3-Clause License. /*! diff --git a/hoomd/mpcd/ReverseNonequilibriumShearFlow.h b/hoomd/mpcd/ReverseNonequilibriumShearFlow.h index f9adab41ae..2a146f4701 100644 --- a/hoomd/mpcd/ReverseNonequilibriumShearFlow.h +++ b/hoomd/mpcd/ReverseNonequilibriumShearFlow.h @@ -1,4 +1,4 @@ -// Copyright (c) 2009-2024 The Regents of the University of Michigan. +// Copyright (c) 2009-2025 The Regents of the University of Michigan. // Part of HOOMD-blue, released under the BSD 3-Clause License. /*! diff --git a/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cc b/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cc index 00ecbc6a88..d4d1d5a24a 100644 --- a/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cc +++ b/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cc @@ -1,4 +1,4 @@ -// Copyright (c) 2009-2024 The Regents of the University of Michigan. +// Copyright (c) 2009-2025 The Regents of the University of Michigan. // Part of HOOMD-blue, released under the BSD 3-Clause License. /*! diff --git a/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cu b/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cu index add25b12f8..230dad1419 100644 --- a/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cu +++ b/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cu @@ -1,4 +1,4 @@ -// Copyright (c) 2009-2024 The Regents of the University of Michigan. +// Copyright (c) 2009-2025 The Regents of the University of Michigan. // Part of HOOMD-blue, released under the BSD 3-Clause License. /*! diff --git a/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cuh b/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cuh index 681a51e514..7742de0a67 100644 --- a/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cuh +++ b/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.cuh @@ -1,4 +1,4 @@ -// Copyright (c) 2009-2024 The Regents of the University of Michigan. +// Copyright (c) 2009-2025 The Regents of the University of Michigan. // Part of HOOMD-blue, released under the BSD 3-Clause License. #ifndef MPCD_REVERSE_NONEQUILIBRIUM_SHEAR_FLOW_GPU_CUH_ diff --git a/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.h b/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.h index 8b7c930fc3..0a9c7bcf27 100644 --- a/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.h +++ b/hoomd/mpcd/ReverseNonequilibriumShearFlowGPU.h @@ -1,4 +1,4 @@ -// Copyright (c) 2009-2024 The Regents of the University of Michigan. +// Copyright (c) 2009-2025 The Regents of the University of Michigan. // Part of HOOMD-blue, released under the BSD 3-Clause License. /*! diff --git a/hoomd/mpcd/ReverseNonequilibriumShearFlowUtilities.h b/hoomd/mpcd/ReverseNonequilibriumShearFlowUtilities.h index 9cfbe04c9b..e88e57f24d 100644 --- a/hoomd/mpcd/ReverseNonequilibriumShearFlowUtilities.h +++ b/hoomd/mpcd/ReverseNonequilibriumShearFlowUtilities.h @@ -1,4 +1,4 @@ -// Copyright (c) 2009-2024 The Regents of the University of Michigan. +// Copyright (c) 2009-2025 The Regents of the University of Michigan. // Part of HOOMD-blue, released under the BSD 3-Clause License. /*! diff --git a/hoomd/mpcd/pytest/test_update.py b/hoomd/mpcd/pytest/test_update.py index 5919ac0c74..9e3e692c6c 100644 --- a/hoomd/mpcd/pytest/test_update.py +++ b/hoomd/mpcd/pytest/test_update.py @@ -1,4 +1,4 @@ -# Copyright (c) 2009-2024 The Regents of the University of Michigan. +# Copyright (c) 2009-2025 The Regents of the University of Michigan. # Part of HOOMD-blue, released under the BSD 3-Clause License. import pytest diff --git a/hoomd/mpcd/update.py b/hoomd/mpcd/update.py index 43e18c36ee..cf9b6bcf19 100644 --- a/hoomd/mpcd/update.py +++ b/hoomd/mpcd/update.py @@ -1,4 +1,4 @@ -# Copyright (c) 2009-2024 The Regents of the University of Michigan. +# Copyright (c) 2009-2025 The Regents of the University of Michigan. # Part of HOOMD-blue, released under the BSD 3-Clause License. r"""MPCD updaters.