Skip to content

Commit

Permalink
Merge pull request #717 from streeve/neighbor_histogram
Browse files Browse the repository at this point in the history 
Create neighbor list histograms
  • Loading branch information
@streeve
streeve authored Dec 7, 2023
2 parents bb5332a + 5d58d8d commit 973e7d0
Show file tree
Hide file tree
Showing 3 changed files with 176 additions and 0 deletions.
51 changes: 51 additions & 0 deletions core/src/Cabana_NeighborList.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,8 @@

#include <Kokkos_Core.hpp>

#include <Cabana_Sort.hpp>

namespace Cabana
{
//---------------------------------------------------------------------------//
Expand Down Expand Up @@ -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 <class ExecutionSpace, class ListType>
Kokkos::View<int* [2], typename ListType::memory_space>
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<int*, typename ListType::memory_space>(
"particle_count", num_particles );

// Extract from neighbor list interface.
auto extract_functor = KOKKOS_LAMBDA( const int p )
{
num_neigh( p ) = NeighborList<ListType>::numNeighbor( list, p );
};
Kokkos::RangePolicy<ExecutionSpace> 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<int* [2], typename ListType::memory_space>(
"particle_count", num_bin, 2 );
auto histogram_functor = KOKKOS_LAMBDA( const int b )
{
int max_neigh = NeighborList<ListType>::maxNeighbor( list );
double bin_width =
static_cast<double>( max_neigh ) / static_cast<double>( 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<int>( ( b + 1 ) * bin_width );
histogram( b, 1 ) = bin_data.binSize( b );
};
Kokkos::RangePolicy<ExecutionSpace> 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
57 changes: 57 additions & 0 deletions core/unit_test/neighbor_unit_test.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<double[3]>;
using AoSoA_t = Cabana::AoSoA<DataTypes, TEST_MEMSPACE>;
AoSoA_t aosoa;

TestNeighborList<typename TEST_EXECSPACE::array_layout, Kokkos::HostSpace>
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<TEST_EXECSPACE> 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
68 changes: 68 additions & 0 deletions core/unit_test/tstNeighborList.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -286,6 +286,7 @@ void testNeighborParallelReduce()
test_data.aosoa, false );
}

//---------------------------------------------------------------------------//
template <class LayoutTag>
void testModifyNeighbors()
{
Expand Down Expand Up @@ -323,6 +324,63 @@ void testModifyNeighbors()
EXPECT_EQ( list_copy.neighbors( p, n ), new_id );
}
}

//---------------------------------------------------------------------------//
template <class LayoutTag>
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<TEST_MEMSPACE, Cabana::FullNeighborTag,
LayoutTag, Cabana::TeamOpTag>;
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
//---------------------------------------------------------------------------//
Expand Down Expand Up @@ -399,6 +457,16 @@ TEST( TEST_CATEGORY, modify_list_test )
#endif
testModifyNeighbors<Cabana::VerletLayout2D>();
}

//---------------------------------------------------------------------------//
TEST( TEST_CATEGORY, neighbor_histogram_test )
{
#ifndef KOKKOS_ENABLE_OPENMPTARGET
testNeighborHistogram<Cabana::VerletLayoutCSR>();
#endif
testNeighborHistogram<Cabana::VerletLayout2D>();
}

//---------------------------------------------------------------------------//

} // end namespace Test

0 comments on commit 973e7d0

Please sign in to comment.