Merge pull request #812 from streeve/reciprocal_neighbors

Guarantee reciprocal neighbors

Author: streeve

core/src/Cabana_VerletList.hpp

@@ -235,16 +235,69 @@ struct VerletListBuilder
         rsqr = neighborhood_radius * neighborhood_radius;
     }
 
-    KOKKOS_INLINE_FUNCTION auto cutoff( [[maybe_unused]] const int p ) const
+    // Check if particle pair i-j is within cutoff, potentially with variable
+    // radii.
+    KOKKOS_INLINE_FUNCTION auto withinCutoff( [[maybe_unused]] const int i,
+                                              const double dist_sqr ) const
     {
         // Square the radius on the fly if using a per-particle field to avoid a
         // deep copy.
         if constexpr ( is_slice<RadiusType>::value ||
                        Kokkos::is_view<RadiusType>::value )
-            return radius( p ) * radius( p );
+            return dist_sqr <= radius( i ) * radius( i );
         // This value is already squared.
         else
-            return rsqr;
+            return dist_sqr <= rsqr;
+    }
+
+    // Check if potential neighbor j is NOT within cutoff, meaning particle i
+    // should add instead for symmetry.
+    KOKKOS_INLINE_FUNCTION auto
+    neighborNotWithinCutoff( [[maybe_unused]] const int j,
+                             [[maybe_unused]] const double dist_sqr ) const
+    {
+        // This neighbor needs to be added if they will not find this particle.
+        if constexpr ( is_slice<RadiusType>::value ||
+                       Kokkos::is_view<RadiusType>::value )
+        {
+            return dist_sqr >= radius( j ) * radius( j );
+        }
+        else
+        {
+            // For a fixed radius, this will never occur.
+            return false;
+        }
+    }
+
+    // Count neighbors, with consideration for self particle i and neighbor j.
+    KOKKOS_INLINE_FUNCTION auto countNeighbor( const int i, const int j,
+                                               const double dist_sqr ) const
+    {
+        int c = 0;
+        // Always add self if within cutoff.
+        if ( withinCutoff( i, dist_sqr ) )
+        {
+            c++;
+            // Add neighbor if they will not find this particle.
+            if ( neighborNotWithinCutoff( j, dist_sqr ) )
+                c++;
+        }
+        return c;
+    }
+
+    // Add neighbors, with consideration for self particle i and neighbor j.
+    KOKKOS_INLINE_FUNCTION void addNeighbor( const int i, const int j,
+                                             const double dist_sqr ) const
+    {
+        // Always add self if within cutoff.
+        if ( withinCutoff( i, dist_sqr ) )
+        {
+            _data.addNeighbor( i, j );
+
+            // Add neighbor if they will not find this particle.
+            if ( neighborNotWithinCutoff( j, dist_sqr ) )
+                _data.addNeighbor( j, i );
+        }
     }
 
     // Neighbor count team operator (only used for CSR lists).
@@ -295,9 +348,11 @@ struct VerletListBuilder
                             {
                                 // See if we should actually check this box for
                                 // neighbors.
-                                if ( linked_cell_list.cellStencil()
-                                         .grid.minDistanceToPoint(
-                                             x_p, y_p, z_p, i, j, k ) <= rsqr )
+                                if ( withinCutoff(
+                                         pid,
+                                         linked_cell_list.cellStencil()
+                                             .grid.minDistanceToPoint(
+                                                 x_p, y_p, z_p, i, j, k ) ) )
                                 {
                                     std::size_t n_offset =
                                         linked_cell_list.binOffset( i, j, k );
@@ -373,8 +428,7 @@ struct VerletListBuilder
             PositionValueType dist_sqr = dx * dx + dy * dy + dz * dz;
 
             // If within the cutoff add to the count.
-            if ( dist_sqr <= cutoff( pid ) )
-                local_count += 1;
+            local_count += countNeighbor( pid, nid, dist_sqr );
         }
     }
 
@@ -512,9 +566,11 @@ struct VerletListBuilder
                             {
                                 // See if we should actually check this box for
                                 // neighbors.
-                                if ( linked_cell_list.cellStencil()
-                                         .grid.minDistanceToPoint(
-                                             x_p, y_p, z_p, i, j, k ) <= rsqr )
+                                if ( withinCutoff(
+                                         pid,
+                                         linked_cell_list.cellStencil()
+                                             .grid.minDistanceToPoint(
+                                                 x_p, y_p, z_p, i, j, k ) ) )
                                 {
                                     // Check the particles in this bin to see if
                                     // they are neighbors.
@@ -583,10 +639,7 @@ struct VerletListBuilder
 
             // If within the cutoff increment the neighbor count and add as a
             // neighbor at that index.
-            if ( dist_sqr <= cutoff( pid ) )
-            {
-                _data.addNeighbor( pid, nid );
-            }
+            addNeighbor( pid, nid, dist_sqr );
         }
     }
 };

core/unit_test/tstNeighborList.hpp

@@ -214,7 +214,7 @@ void testNonUniformRadius()
     // Create the AoSoA and fill custom particle details.
     std::size_t particle_x = 2;
     // Purposely choose radius to reach all particles.
-    double large_radius = 4.35;
+    double large_radius = 4.05;
     // Purposely choose radius to reach nearest neighbors only.
     double small_radius = 3.32;
     // Create the AoSoA and fill with particles
@@ -249,7 +249,7 @@ void testNonUniformRadius()
         if ( p == 0 || p == test_data.num_particle - 1 )
             EXPECT_EQ( list_copy.counts( p ), 6 );
         else
-            EXPECT_EQ( list_copy.counts( p ), 3 );
+            EXPECT_EQ( list_copy.counts( p ), 4 );
     }
 }