diff --git a/core/src/Cabana_NeighborList.hpp b/core/src/Cabana_NeighborList.hpp index ca05faf95..934f04d4e 100644 --- a/core/src/Cabana_NeighborList.hpp +++ b/core/src/Cabana_NeighborList.hpp @@ -18,6 +18,8 @@ #include +#include + namespace Cabana { //---------------------------------------------------------------------------// @@ -114,6 +116,55 @@ maxNeighbor( const ListType& list, const std::size_t num_particles ) } } // namespace Impl +/*! + \brief Construct a histogram of neighbors per particle. + \param exec_space Kokkos execution space. + \param num_particles Number of particles. + \param list Neighbor list with valid NeighborList interface. + \param num_bin Number of bins for the histogram. + \return Neighbor list histogram View. +*/ +template +Kokkos::View +neighborHistogram( ExecutionSpace exec_space, const std::size_t num_particles, + const ListType& list, const int num_bin ) +{ + // Allocate View of neighbors per particle + auto num_neigh = Kokkos::View( + "particle_count", num_particles ); + + // Extract from neighbor list interface. + auto extract_functor = KOKKOS_LAMBDA( const int p ) + { + num_neigh( p ) = NeighborList::numNeighbor( list, p ); + }; + Kokkos::RangePolicy particle_policy( exec_space, 0, + num_particles ); + Kokkos::parallel_for( particle_policy, extract_functor ); + Kokkos::fence(); + + auto bin_data = Cabana::binByKey( num_neigh, num_bin ); + + auto histogram = Kokkos::View( + "particle_count", num_bin, 2 ); + auto histogram_functor = KOKKOS_LAMBDA( const int b ) + { + int max_neigh = NeighborList::maxNeighbor( list ); + double bin_width = + static_cast( max_neigh ) / static_cast( num_bin ); + if ( num_bin > max_neigh ) + bin_width = 1; + // Wait to cast back to int to get the actual bin edge. + histogram( b, 0 ) = static_cast( ( b + 1 ) * bin_width ); + histogram( b, 1 ) = bin_data.binSize( b ); + }; + Kokkos::RangePolicy bin_policy( exec_space, 0, num_bin ); + Kokkos::parallel_for( bin_policy, histogram_functor ); + Kokkos::fence(); + + return histogram; +} + } // end namespace Cabana #endif // end CABANA_NEIGHBORLIST_HPP diff --git a/core/unit_test/neighbor_unit_test.hpp b/core/unit_test/neighbor_unit_test.hpp index 3633db808..162e3f6ff 100644 --- a/core/unit_test/neighbor_unit_test.hpp +++ b/core/unit_test/neighbor_unit_test.hpp @@ -1001,5 +1001,62 @@ struct NeighborListTestData } }; +//---------------------------------------------------------------------------// +// Default ordered test settings. +struct NeighborListTestDataOrdered +{ + int num_particle; + int num_ignore = 100; + double test_radius; + double box_min = 0.0; + double box_max = 5.0; + + double cell_size_ratio = 0.5; + double grid_min[3] = { box_min, box_min, box_min }; + double grid_max[3] = { box_max, box_max, box_max }; + + using DataTypes = Cabana::MemberTypes; + using AoSoA_t = Cabana::AoSoA; + AoSoA_t aosoa; + + TestNeighborList + N2_list_copy; + + NeighborListTestDataOrdered( const int particle_x, const int m = 3 ) + { + num_particle = particle_x * particle_x * particle_x; + double dx = ( grid_max[0] - grid_min[0] ) / particle_x; + // Use a fixed ratio of cutoff / spacing (and include a floating point + // tolerance). + test_radius = m * dx + 1e-7; + + // Create the AoSoA and fill with ordered particle positions. + aosoa = AoSoA_t( "ordered", num_particle ); + createParticles( particle_x, dx ); + + // Create a full N^2 neighbor list to check against. + auto positions = Cabana::slice<0>( aosoa ); + auto N2_list = computeFullNeighborList( positions, test_radius ); + N2_list_copy = createTestListHostCopy( N2_list ); + } + + void createParticles( const int particle_x, const double dx ) + { + auto positions = Cabana::slice<0>( aosoa ); + + Kokkos::RangePolicy policy( 0, positions.size() ); + Kokkos::parallel_for( + "ordered_particles", policy, KOKKOS_LAMBDA( int pid ) { + int i = pid / ( particle_x * particle_x ); + int j = ( pid / particle_x ) % particle_x; + int k = pid % particle_x; + positions( pid, 0 ) = dx / 2 + dx * i; + positions( pid, 1 ) = dx / 2 + dx * j; + positions( pid, 2 ) = dx / 2 + dx * k; + } ); + Kokkos::fence(); + } +}; + //---------------------------------------------------------------------------// } // end namespace Test diff --git a/core/unit_test/tstNeighborList.hpp b/core/unit_test/tstNeighborList.hpp index 1bb28fd8c..25301bb89 100644 --- a/core/unit_test/tstNeighborList.hpp +++ b/core/unit_test/tstNeighborList.hpp @@ -286,6 +286,7 @@ void testNeighborParallelReduce() test_data.aosoa, false ); } +//---------------------------------------------------------------------------// template void testModifyNeighbors() { @@ -323,6 +324,63 @@ void testModifyNeighbors() EXPECT_EQ( list_copy.neighbors( p, n ), new_id ); } } + +//---------------------------------------------------------------------------// +template +void testNeighborHistogram() +{ + int particle_x = 10; + // Create the AoSoA and fill with particles + NeighborListTestDataOrdered test_data( particle_x ); + auto position = Cabana::slice<0>( test_data.aosoa ); + + // Create the neighbor list. + using ListType = Cabana::VerletList; + ListType nlist( position, 0, position.size(), test_data.test_radius, + test_data.cell_size_ratio, test_data.grid_min, + test_data.grid_max ); + + // Create the neighbor histogram + { + int num_bin = 10; + auto nhist = neighborHistogram( TEST_EXECSPACE{}, position.size(), + nlist, num_bin ); + auto host_histogram = + Kokkos::create_mirror_view_and_copy( Kokkos::HostSpace(), nhist ); + + // 122 is the number of neighbors expected for a full spherical shell + // with characteristic ratio (cutoff / dx) = 3 + double bin_max[10] = { 12, 24, 36, 48, 61, 73, 85, 97, 109, 122 }; + // This is the direct copied output (zeros due to fixed spacing). + double bin_count[10] = { 32, 72, 24, 152, 120, 168, 0, 216, 0, 152 }; + + for ( int i = 0; i < num_bin; ++i ) + { + EXPECT_EQ( bin_max[i], host_histogram( i, 0 ) ); + EXPECT_EQ( bin_count[i], host_histogram( i, 1 ) ); + } + } + + // Create another histogram with fewer bins. + { + int num_bin = 5; + auto nhist = neighborHistogram( TEST_EXECSPACE{}, position.size(), + nlist, num_bin ); + auto host_histogram = + Kokkos::create_mirror_view_and_copy( Kokkos::HostSpace(), nhist ); + + double bin_max[5] = { 24, 48, 73, 97, 122 }; + // This is the direct copied output. + double bin_count[5] = { 104, 176, 288, 216, 152 }; + for ( int i = 0; i < num_bin; ++i ) + { + EXPECT_EQ( bin_max[i], host_histogram( i, 0 ) ); + EXPECT_EQ( bin_count[i], host_histogram( i, 1 ) ); + } + } +} + //---------------------------------------------------------------------------// // TESTS //---------------------------------------------------------------------------// @@ -399,6 +457,16 @@ TEST( TEST_CATEGORY, modify_list_test ) #endif testModifyNeighbors(); } + +//---------------------------------------------------------------------------// +TEST( TEST_CATEGORY, neighbor_histogram_test ) +{ +#ifndef KOKKOS_ENABLE_OPENMPTARGET + testNeighborHistogram(); +#endif + testNeighborHistogram(); +} + //---------------------------------------------------------------------------// } // end namespace Test