Lees-Edwards boundary conditions for shear flow

hoomd/BoxResizeUpdater.cc

@@ -96,6 +96,9 @@ void BoxResizeUpdater::scaleAndWrapParticles(const BoxDim& cur_box, const BoxDim
 
     // ensure that the particles are still in their
     // local boxes by wrapping them if they are not
+    ArrayHandle<Scalar4> h_vel(m_pdata->getVelocities(),
+                               access_location::host,
+                               access_mode::readwrite);
     ArrayHandle<int3> h_image(m_pdata->getImages(), access_location::host, access_mode::readwrite);
 
     const BoxDim& local_box = m_pdata->getBox();
@@ -104,7 +107,7 @@ void BoxResizeUpdater::scaleAndWrapParticles(const BoxDim& cur_box, const BoxDim
         {
         // need to update the image if we move particles from one side
         // of the box to the other
-        local_box.wrap(h_pos.data[i], h_image.data[i]);
+        local_box.wrap(h_pos.data[i], h_vel.data[i], h_image.data[i]);
         }
     }
 

hoomd/BoxResizeUpdaterGPU.cc

@@ -42,6 +42,10 @@ void BoxResizeUpdaterGPU::scaleAndWrapParticles(const BoxDim& cur_box, const Box
                                access_location::device,
                                access_mode::readwrite);
 
+    ArrayHandle<Scalar4> d_vel(m_pdata->getVelocities(),
+                               access_location::device,
+                               access_mode::readwrite);
+
     ArrayHandle<int3> d_image(m_pdata->getImages(),
                               access_location::device,
                               access_mode::readwrite);
@@ -62,6 +66,7 @@ void BoxResizeUpdaterGPU::scaleAndWrapParticles(const BoxDim& cur_box, const Box
     m_tuner_wrap->begin();
     kernel::gpu_box_resize_wrap(m_pdata->getN(),
                                 d_pos.data,
+                                d_vel.data,
                                 d_image.data,
                                 new_box,
                                 m_tuner_wrap->getParam()[0]);

hoomd/BoxResizeUpdaterGPU.cu

@@ -31,14 +31,17 @@ __global__ void gpu_box_resize_scale_kernel(Scalar4* d_pos,
         }
     }
 
-__global__ void
-gpu_box_resize_wrap_kernel(unsigned int N, Scalar4* d_pos, int3* d_image, const BoxDim new_box)
+__global__ void gpu_box_resize_wrap_kernel(unsigned int N,
+                                           Scalar4* d_pos,
+                                           Scalar4* d_vel,
+                                           int3* d_image,
+                                           const BoxDim new_box)
     {
     int idx = blockIdx.x * blockDim.x + threadIdx.x;
 
     if (idx < N)
         {
-        new_box.wrap(d_pos[idx], d_image[idx]);
+        new_box.wrap(d_pos[idx], d_vel[idx], d_image[idx]);
         }
     }
 
@@ -74,6 +77,7 @@ hipError_t gpu_box_resize_scale(Scalar4* d_pos,
 
 hipError_t gpu_box_resize_wrap(const unsigned int N,
                                Scalar4* d_pos,
+                               Scalar4* d_vel,
                                int3* d_image,
                                const BoxDim& new_box,
                                unsigned int block_size)
@@ -94,6 +98,7 @@ hipError_t gpu_box_resize_wrap(const unsigned int N,
                        0,
                        N,
                        d_pos,
+                       d_vel,
                        d_image,
                        new_box);
 

hoomd/BoxResizeUpdaterGPU.cuh

@@ -21,6 +21,7 @@ hipError_t gpu_box_resize_scale(Scalar4* d_pos,
 
 hipError_t gpu_box_resize_wrap(const unsigned int N,
                                Scalar4* d_pos,
+                               Scalar4* d_vel,
                                int3* d_image,
                                const BoxDim& new_box,
                                unsigned int block_size);

hoomd/Communicator.cc

@@ -1719,9 +1719,10 @@ void Communicator::migrateParticles()
             {
             detail::pdata_element& p = m_recvbuf[idx];
             Scalar4& postype = p.pos;
+            Scalar4& vel = p.vel;
             int3& image = p.image;
 
-            shifted_box.wrap(postype, image);
+            shifted_box.wrap(postype, vel, image);
             }
 
         // remove particles that were sent and fill particle data with received particles
@@ -2389,11 +2390,14 @@ void Communicator::exchangeGhosts()
             }
 
         // wrap particle positions
-        if (flags[comm_flag::position])
+        if (flags[comm_flag::position] || flags[comm_flag::velocity])
             {
             ArrayHandle<Scalar4> h_pos(m_pdata->getPositions(),
                                        access_location::host,
                                        access_mode::readwrite);
+            ArrayHandle<Scalar4> h_vel(m_pdata->getVelocities(),
+                                       access_location::host,
+                                       access_mode::readwrite);
             ArrayHandle<int3> h_image(m_pdata->getImages(),
                                       access_location::host,
                                       access_mode::readwrite);
@@ -2403,10 +2407,18 @@ void Communicator::exchangeGhosts()
             for (unsigned int idx = start_idx; idx < start_idx + m_num_recv_ghosts[dir]; idx++)
                 {
                 Scalar4& pos = h_pos.data[idx];
+                int3& img = h_image.data[idx];
 
                 // wrap particles received across a global boundary
-                int3& img = h_image.data[idx];
-                shifted_box.wrap(pos, img);
+                if (flags[comm_flag::velocity])
+                    {
+                    Scalar4& vel = h_vel.data[idx];
+                    shifted_box.wrap(pos, vel, img);
+                    }
+                else
+                    {
+                    shifted_box.wrap(pos, img);
+                    }
                 }
             }
 
@@ -2961,20 +2973,31 @@ void Communicator::beginUpdateGhosts(uint64_t timestep)
             MPI_Waitall(2, &m_reqs.front(), &m_stats.front());
             }
         // wrap particle positions (only if copying positions)
-        if (flags[comm_flag::position])
+        if (flags[comm_flag::position] || flags[comm_flag::velocity])
             {
             ArrayHandle<Scalar4> h_pos(m_pdata->getPositions(),
                                        access_location::host,
                                        access_mode::readwrite);
+            ArrayHandle<Scalar4> h_vel(m_pdata->getVelocities(),
+                                       access_location::host,
+                                       access_mode::readwrite);
 
             const BoxDim shifted_box = getShiftedBox();
             for (unsigned int idx = start_idx; idx < start_idx + m_num_recv_ghosts[dir]; idx++)
                 {
                 Scalar4& pos = h_pos.data[idx];
+                int3 img = make_int3(0, 0, 0);
 
                 // wrap particles received across a global boundary
-                int3 img = make_int3(0, 0, 0);
-                shifted_box.wrap(pos, img);
+                if (flags[comm_flag::velocity])
+                    {
+                    Scalar4& vel = h_vel.data[idx];
+                    shifted_box.wrap(pos, vel, img);
+                    }
+                else
+                    {
+                    shifted_box.wrap(pos, img);
+                    }
                 }
             }
 

hoomd/CommunicatorGPU.cu

@@ -260,7 +260,7 @@ __global__ void gpu_wrap_particles_kernel(const unsigned int n_recv,
         return;
 
     detail::pdata_element p = d_recv[idx];
-    box.wrap(p.pos, p.image);
+    box.wrap(p.pos, p.vel, p.image);
     d_recv[idx] = p;
     }
 
@@ -789,8 +789,10 @@ gpu_pack_kernel(unsigned int n_out, const uint2* d_ghost_idx_adj, const T* in, T
 __global__ void gpu_pack_wrap_kernel(unsigned int n_out,
                                      const uint2* d_ghost_idx_adj,
                                      const Scalar4* d_postype,
+                                     const Scalar4* d_vel,
                                      const int3* d_img,
                                      Scalar4* out_pos,
+                                     Scalar4* out_vel,
                                      int3* out_img,
                                      Index3D di,
                                      uint3 my_pos,
@@ -860,14 +862,24 @@ __global__ void gpu_pack_wrap_kernel(unsigned int n_out,
         }
 
     box.setPeriodic(periodic);
+    Scalar4 postype = d_postype[idx];
     int3 img = make_int3(0, 0, 0);
-    if (d_img)
+    if (out_img)
+        {
         img = d_img[idx];
-    Scalar4 postype = d_postype[idx];
-    box.wrap(postype, img, wrap);
+        }
+    if (out_vel)
+        {
+        Scalar4 vel = d_vel[idx];
+        box.wrap(postype, vel, img, wrap);
+        out_vel[buf_idx] = vel;
+        }
+    else
+        {
+        box.wrap(postype, img, wrap);
+        }
 
     out_pos[buf_idx] = postype;
-
     if (out_img)
         {
         out_img[buf_idx] = img;
@@ -928,10 +940,19 @@ void gpu_exchange_ghosts_pack(unsigned int n_out,
                            d_tag,
                            d_tag_sendbuf);
         }
-    if (send_pos)
+    /*
+     * combined packing pathway for positions , velocities and images
+     * call gpu_pack_wrap_kernel whenever any of send_pos, send_vel and send_image is true
+     * for buffers that are not being sent, pass a nullptr to the respective output
+     */
+    if (send_pos || send_vel || send_image)
         {
         assert(d_pos);
         assert(d_pos_sendbuf);
+        if (send_vel)
+            {
+            assert(d_vel);
+            }
         if (send_image)
             {
             assert(d_img);
@@ -944,27 +965,15 @@ void gpu_exchange_ghosts_pack(unsigned int n_out,
                            n_out,
                            d_ghost_idx_adj,
                            d_pos,
+                           d_vel,
                            d_img,
                            d_pos_sendbuf,
+                           send_vel ? d_vel_sendbuf : 0,
                            send_image ? d_img_sendbuf : 0,
                            di,
                            my_pos,
                            box);
         }
-    if (send_vel)
-        {
-        assert(d_vel);
-        assert(d_vel_sendbuf);
-        hipLaunchKernelGGL(gpu_pack_kernel,
-                           dim3(n_blocks),
-                           dim3(block_size),
-                           0,
-                           0,
-                           n_out,
-                           d_ghost_idx_adj,
-                           d_vel,
-                           d_vel_sendbuf);
-        }
     if (send_charge)
         {
         assert(d_charge);

hoomd/GSDDumpWriter.cc

@@ -1113,6 +1113,9 @@ void GSDDumpWriter::populateLocalFrame(GSDDumpWriter::GSDFrame& frame, uint64_t
         ArrayHandle<Scalar4> h_postype(m_pdata->getPositions(),
                                        access_location::host,
                                        access_mode::read);
+        ArrayHandle<Scalar4> h_veltype(m_pdata->getVelocities(),
+                                       access_location::host,
+                                       access_mode::read);
         ArrayHandle<int3> h_image(m_pdata->getImages(), access_location::host, access_mode::read);
 
         if (m_dynamic[gsd_flag::particles_position] || m_nframes == 0)
@@ -1132,6 +1135,7 @@ void GSDDumpWriter::populateLocalFrame(GSDDumpWriter::GSDFrame& frame, uint64_t
             {
             vec3<Scalar> position
                 = vec3<Scalar>(h_postype.data[index]) - vec3<Scalar>(m_pdata->getOrigin());
+            vec3<Scalar> velocity = vec3<Scalar>(h_veltype.data[index]);
             unsigned int type = __scalar_as_int(h_postype.data[index].w);
             int3 image = make_int3(0, 0, 0);
 
@@ -1140,7 +1144,7 @@ void GSDDumpWriter::populateLocalFrame(GSDDumpWriter::GSDFrame& frame, uint64_t
                 image = h_image.data[index];
                 }
 
-            frame.global_box.wrap(position, image);
+            frame.global_box.wrap(position, velocity, image);
 
             if (m_dynamic[gsd_flag::particles_position] || m_nframes == 0)
                 {