diff --git a/core/src/Cabana_N2NeighborList.hpp b/core/src/Cabana_N2NeighborList.hpp new file mode 100644 index 000000000..9b5451e66 --- /dev/null +++ b/core/src/Cabana_N2NeighborList.hpp @@ -0,0 +1,602 @@ +/**************************************************************************** + * Copyright (c) 2018-2022 by the Cabana authors * + * All rights reserved. * + * * + * This file is part of the Cabana library. Cabana is distributed under a * + * BSD 3-clause license. For the licensing terms see the LICENSE file in * + * the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +/*! + \file Cabana_N2NeighborList.hpp + \brief Neighbor list without grid acceleration +*/ +#ifndef CABANA_N2LIST_HPP +#define CABANA_N2LIST_HPP + +#include +#include + +#include + +#include + +namespace Cabana +{ + +namespace Impl +{ +//! \cond Impl + +struct N2ListTag +{ +}; + +//---------------------------------------------------------------------------// +// Neighbor List Builder +//---------------------------------------------------------------------------// +template +struct N2NeighborListBuilder +{ + // Types. + using device = DeviceType; + using PositionValueType = typename PositionSlice::value_type; + using RandomAccessPositionSlice = + typename PositionSlice::random_access_slice; + using memory_space = typename device::memory_space; + using execution_space = typename device::execution_space; + + // List data. + NeighborListData _data; + + // Neighbor cutoff. + PositionValueType rsqr; + + // Positions. + RandomAccessPositionSlice position; + std::size_t pid_begin, pid_end; + + // Check to count or refill. + bool refill; + bool count; + + // Maximum neighbors per particle + std::size_t max_n; + + // Constructor. + N2NeighborListBuilder( PositionSlice slice, const std::size_t begin, + const std::size_t end, + const PositionValueType neighborhood_radius, + const std::size_t max_neigh ) + : pid_begin( begin ) + , pid_end( end ) + , max_n( max_neigh ) + { + count = true; + refill = false; + + // Create the count view. + _data.counts = + Kokkos::View( "num_neighbors", slice.size() ); + + // Make a guess for the number of neighbors per particle for 2D + // lists. + initCounts( LayoutTag() ); + + // Get the positions with random access read-only memory. + position = slice; + + // We will use the square of the distance for neighbor + // determination. + rsqr = neighborhood_radius * neighborhood_radius; + } + + // Neighbor count team operator (only used for CSR lists). + struct CountNeighborsTag + { + }; + KOKKOS_INLINE_FUNCTION + void operator()( const CountNeighborsTag&, const int pid ) const + { + if ( ( pid >= pid_begin ) && ( pid < pid_end ) ) + { + // Cache the particle coordinates. + double x_p = position( pid, 0 ); + double y_p = position( pid, 1 ); + double z_p = position( pid, 2 ); + + // Check to see if the particles are neighbors. + int count = 0; + // Note that we loop over all particles as potential neighbors + // of particles in the given range. + for ( std::size_t nid = 0; nid < position.size(); nid++ ) + // Check to see if the particles are neighbors. + { + neighbor_kernel( pid, x_p, y_p, z_p, nid, count ); + } + _data.counts( pid ) = count; + } + } + + using CountNeighborsPolicy = + Kokkos::TeamPolicy, + Kokkos::Schedule>; + KOKKOS_INLINE_FUNCTION + void + operator()( const CountNeighborsTag&, + const typename CountNeighborsPolicy::member_type& team ) const + { + // The league rank of the team is the current particle. + std::size_t pid = team.league_rank() + pid_begin; + if ( ( pid >= pid_begin ) && ( pid < pid_end ) ) + { + // Cache the particle coordinates. + double x_p = position( pid, 0 ); + double y_p = position( pid, 1 ); + double z_p = position( pid, 2 ); + + // Check to see if the particles are neighbors. + int count = 0; + // Note that we loop over all particles as potential neighbors + // of particles in the given range. + Kokkos::parallel_reduce( + Kokkos::TeamThreadRange( team, position.size() ), + [&]( const int nid, int& local_count ) + { + // Check to see if the particles are neighbors. + neighbor_kernel( pid, x_p, y_p, z_p, nid, local_count ); + }, + count ); + _data.counts( pid ) = count; + } + } + + // Neighbor count kernel + KOKKOS_INLINE_FUNCTION + void neighbor_kernel( const int pid, const double x_p, const double y_p, + const double z_p, const int nid, + int& local_count ) const + { + // Cache the candidate neighbor particle coordinates. + double x_n = position( nid, 0 ); + double y_n = position( nid, 1 ); + double z_n = position( nid, 2 ); + + // If this could be a valid neighbor, continue. + if ( NeighborDiscriminator::isValid( + pid, x_p, y_p, z_p, nid, x_n, y_n, z_n ) ) + { + // Calculate the distance between the particle and its candidate + // neighbor. + PositionValueType dx = x_p - x_n; + PositionValueType dy = y_p - y_n; + PositionValueType dz = z_p - z_n; + PositionValueType dist_sqr = dx * dx + dy * dy + dz * dz; + + // If within the cutoff add to the count. + if ( dist_sqr <= rsqr ) + local_count += 1; + } + } + + // Process the CSR counts by computing offsets and allocating the + // neighbor list. + template + struct OffsetScanOp + { + using kokkos_mem_space = KokkosMemorySpace; + Kokkos::View counts; + Kokkos::View offsets; + KOKKOS_INLINE_FUNCTION + void operator()( const int i, int& update, const bool final_pass ) const + { + if ( final_pass ) + offsets( i ) = update; + update += counts( i ); + } + }; + + void initCounts( NeighborLayoutCSR ) {} + + void initCounts( NeighborLayout2D ) + { + if ( max_n > 0 ) + { + count = false; + + _data.neighbors = Kokkos::View( + Kokkos::ViewAllocateWithoutInitializing( "neighbors" ), + _data.counts.size(), max_n ); + } + } + + void processCounts( NeighborLayoutCSR ) + { + // Allocate offsets. + _data.offsets = Kokkos::View( + Kokkos::ViewAllocateWithoutInitializing( "neighbor_offsets" ), + _data.counts.size() ); + + // Calculate offsets from counts and the total number of counts. + OffsetScanOp offset_op; + offset_op.counts = _data.counts; + offset_op.offsets = _data.offsets; + int total_num_neighbor; + Kokkos::RangePolicy range_policy( + 0, _data.counts.extent( 0 ) ); + Kokkos::parallel_scan( "Cabana::NeighborListBuilder::offset_scan", + range_policy, offset_op, total_num_neighbor ); + Kokkos::fence(); + + // Allocate the neighbor list. + _data.neighbors = Kokkos::View( + Kokkos::ViewAllocateWithoutInitializing( "neighbors" ), + total_num_neighbor ); + + // Reset the counts. We count again when we fill. + Kokkos::deep_copy( _data.counts, 0 ); + } + + // Process 2D counts by computing the maximum number of neighbors and + // reallocating the 2D data structure if needed. + void processCounts( NeighborLayout2D ) + { + // Calculate the maximum number of neighbors. + auto counts = _data.counts; + int max_num_neighbor = 0; + Kokkos::Max max_reduce( max_num_neighbor ); + Kokkos::parallel_reduce( + "Cabana::NeighborListBuilder::reduce_max", + Kokkos::RangePolicy( 0, _data.counts.size() ), + KOKKOS_LAMBDA( const int i, int& value ) { + if ( counts( i ) > value ) + value = counts( i ); + }, + max_reduce ); + Kokkos::fence(); + + // Reallocate the neighbor list if previous size is exceeded. + if ( count or ( std::size_t ) + max_num_neighbor > _data.neighbors.extent( 1 ) ) + { + refill = true; + Kokkos::deep_copy( _data.counts, 0 ); + _data.neighbors = Kokkos::View( + Kokkos::ViewAllocateWithoutInitializing( "neighbors" ), + _data.counts.size(), max_num_neighbor ); + } + } + + // Compatibility wrapper for old tags. + void processCounts( VerletLayoutCSR ) + { + processCounts( NeighborLayoutCSR{} ); + } + void processCounts( VerletLayout2D ) + { + processCounts( NeighborLayout2D{} ); + } + void initCounts( VerletLayoutCSR ) { initCounts( NeighborLayoutCSR{} ); } + void initCounts( VerletLayout2D ) { initCounts( NeighborLayout2D{} ); } + + // Neighbor count team operator. + struct FillNeighborsTag + { + }; + KOKKOS_INLINE_FUNCTION + void operator()( const FillNeighborsTag&, const int pid ) const + { + if ( ( pid >= pid_begin ) && ( pid < pid_end ) ) + { + // Cache the particle coordinates. + double x_p = position( pid, 0 ); + double y_p = position( pid, 1 ); + double z_p = position( pid, 2 ); + + // Note that we loop over all particles as potential neighbors + // of particles in the given range. + for ( std::size_t nid = 0; nid < position.size(); nid++ ) + // Check to see if the particles are neighbors. + neighbor_kernel( pid, x_p, y_p, z_p, nid ); + } + }; + + using FillNeighborsPolicy = + Kokkos::TeamPolicy, + Kokkos::Schedule>; + KOKKOS_INLINE_FUNCTION + void + operator()( const FillNeighborsTag&, + const typename FillNeighborsPolicy::member_type& team ) const + { + // The league rank of the team is the current particle. + std::size_t pid = team.league_rank() + pid_begin; + if ( ( pid >= pid_begin ) && ( pid < pid_end ) ) + { + // Cache the particle coordinates. + double x_p = position( pid, 0 ); + double y_p = position( pid, 1 ); + double z_p = position( pid, 2 ); + + // Note that we loop over all particles as potential neighbors + // of particles in the given range. + Kokkos::parallel_for( + Kokkos::TeamThreadRange( team, position.size() ), + [&]( const int nid ) + { + // Check to see if the particles are neighbors. + neighbor_kernel( pid, x_p, y_p, z_p, nid ); + } ); + } + }; + + // Neighbor fill kernel. + KOKKOS_INLINE_FUNCTION + void neighbor_kernel( const int pid, const double x_p, const double y_p, + const double z_p, const int nid ) const + { + // Cache the candidate neighbor particle coordinates. + double x_n = position( nid, 0 ); + double y_n = position( nid, 1 ); + double z_n = position( nid, 2 ); + + // If this could be a valid neighbor, continue. + if ( NeighborDiscriminator::isValid( + pid, x_p, y_p, z_p, nid, x_n, y_n, z_n ) ) + { + // Calculate the distance between the particle and its candidate + // neighbor. + PositionValueType dx = x_p - x_n; + PositionValueType dy = y_p - y_n; + PositionValueType dz = z_p - z_n; + PositionValueType dist_sqr = dx * dx + dy * dy + dz * dz; + + // If within the cutoff increment the neighbor count and add as + // a neighbor at that index. + if ( dist_sqr <= rsqr ) + { + _data.addNeighbor( pid, nid ); + } + } + } +}; + +//---------------------------------------------------------------------------// + +//! \endcond +} // end namespace Impl + +//---------------------------------------------------------------------------// +/*! + \brief Neighbor list implementation based on interaction distance. + + \tparam MemorySpace The Kokkos memory space for storing the neighbor list. + + \tparam AlgorithmTag Tag indicating whether to build a full or half + neighbor list. + + \tparam LayoutTag Tag indicating whether to use a CSR or 2D data layout. + + \tparam BuildTag Tag indicating whether to use hierarchical team or team + vector parallelism when building neighbor lists. +*/ +template +class N2NeighborList +{ + public: + static_assert( Kokkos::is_memory_space::value, "" ); + + //! Kokkos memory space in which the neighbor list data resides. + using memory_space = MemorySpace; + + //! Kokkos default execution space for this memory space. + using execution_space = typename memory_space::execution_space; + + //! Neighbor list data. + NeighborListData _data; + + /*! + \brief Default constructor. + */ + N2NeighborList() {} + + /*! + \brief N2NeighborList constructor. Given a list of particle positions + and a neighborhood radius calculate the neighbor list. + + \param x The slice containing the particle positions + + \param begin The beginning particle index to compute neighbors for. + + \param end The end particle index to compute neighbors for. + + \param neighborhood_radius The radius of the neighborhood. Particles + within this radius are considered neighbors. + + \param max_neigh Optional maximum number of neighbors per particle to + pre-allocate the neighbor list. Potentially avoids recounting with 2D + layout only. + + Particles outside of the neighborhood radius will not be considered + neighbors. Only compute the neighbors of those that are within the + given range. All particles are candidates for being a neighbor, + regardless of whether or not they are in the range. + */ + template + N2NeighborList( + PositionSlice x, const std::size_t begin, const std::size_t end, + const typename PositionSlice::value_type neighborhood_radius, + const std::size_t max_neigh = 0, + typename std::enable_if<( is_slice::value ), + int>::type* = 0 ) + { + build( x, begin, end, neighborhood_radius, max_neigh ); + } + + /*! + \brief Given a list of particle positions and a neighborhood radius + calculate the neighbor list. + */ + template + void build( PositionSlice x, const std::size_t begin, const std::size_t end, + const typename PositionSlice::value_type neighborhood_radius, + const std::size_t max_neigh = 0 ) + { + // Use the default execution space. + build( execution_space{}, x, begin, end, neighborhood_radius, + max_neigh ); + } + /*! + \brief Given a list of particle positions and a neighborhood radius + calculate the neighbor list. + */ + template + void build( ExecutionSpace, PositionSlice x, const std::size_t begin, + const std::size_t end, + const typename PositionSlice::value_type neighborhood_radius, + const std::size_t max_neigh = 0 ) + { + static_assert( is_accessible_from{}, "" ); + + assert( end >= begin ); + assert( end <= x.size() ); + + using device_type = Kokkos::Device; + + // Create a builder functor. + using builder_type = + Impl::N2NeighborListBuilder; + builder_type builder( x, begin, end, neighborhood_radius, max_neigh ); + + // For each particle in the range check for neighbor particles. For CSR + // lists, we count, then fill neighbors. For 2D lists, we count and fill + // at the same time, unless the array size is exceeded, at which point + // only counting is continued to reallocate and refill. + typename builder_type::FillNeighborsPolicy fill_policy( + end - begin, Kokkos::AUTO, 4 ); + if ( builder.count ) + { + typename builder_type::CountNeighborsPolicy count_policy( + end - begin, Kokkos::AUTO, 4 ); + Kokkos::parallel_for( "Cabana::N2NeighborList::count_neighbors", + count_policy, builder ); + } + else + { + builder.processCounts( LayoutTag() ); + Kokkos::parallel_for( "Cabana::N2NeighborList::fill_neighbors", + fill_policy, builder ); + } + Kokkos::fence(); + + // Process the counts by computing offsets and allocating the + // neighbor list, if needed. + builder.processCounts( LayoutTag() ); + + // For each particle in the range fill (or refill) its part of the + // neighbor list. + if ( builder.count or builder.refill ) + { + Kokkos::parallel_for( "Cabana::N2NeighborList::fill_neighbors", + fill_policy, builder ); + Kokkos::fence(); + } + + // Get the data from the builder. + _data = builder._data; + } +}; + +//---------------------------------------------------------------------------// +// Neighbor list interface implementation. +//---------------------------------------------------------------------------// +//! CSR N2 NeighborList interface. +template +class NeighborList< + N2NeighborList> +{ + public: + //! Kokkos memory space. + using memory_space = MemorySpace; + //! Neighbor list type. + using list_type = + N2NeighborList; + + //! Get the total number of neighbors (maximum size of CSR list). + KOKKOS_INLINE_FUNCTION + static std::size_t maxNeighbor( const list_type& list ) + { + return list._data.neighbors.extent( 0 ); + } + + //! Get the number of neighbors for a given particle index. + KOKKOS_INLINE_FUNCTION + static std::size_t numNeighbor( const list_type& list, + const std::size_t particle_index ) + { + return list._data.counts( particle_index ); + } + + //! Get the id for a neighbor for a given particle index and the index + //! of the neighbor relative to the particle. + KOKKOS_INLINE_FUNCTION + static std::size_t getNeighbor( const list_type& list, + const std::size_t particle_index, + const std::size_t neighbor_index ) + { + return list._data.neighbors( list._data.offsets( particle_index ) + + neighbor_index ); + } +}; + +//---------------------------------------------------------------------------// +//! 2D N2 NeighborList interface. +template +class NeighborList< + N2NeighborList> +{ + public: + //! Kokkos memory space. + using memory_space = MemorySpace; + //! Neighbor list type. + using list_type = + N2NeighborList; + + //! Get the maximum number of neighbors per particle. + KOKKOS_INLINE_FUNCTION + static std::size_t maxNeighbor( const list_type& list ) + { + return list._data.neighbors.extent( 1 ); + } + + //! Get the number of neighbors for a given particle index. + KOKKOS_INLINE_FUNCTION + static std::size_t numNeighbor( const list_type& list, + const std::size_t particle_index ) + { + return list._data.counts( particle_index ); + } + + //! Get the id for a neighbor for a given particle index and the index + //! of the neighbor relative to the particle. + KOKKOS_INLINE_FUNCTION + static std::size_t getNeighbor( const list_type& list, + const std::size_t particle_index, + const std::size_t neighbor_index ) + { + return list._data.neighbors( particle_index, neighbor_index ); + } +}; + +//---------------------------------------------------------------------------// + +} // end namespace Cabana + +#endif // end CABANA_N2LIST_HPP diff --git a/core/src/Cabana_VerletList.hpp b/core/src/Cabana_VerletList.hpp index 643d542cc..dde132eb1 100644 --- a/core/src/Cabana_VerletList.hpp +++ b/core/src/Cabana_VerletList.hpp @@ -17,6 +17,7 @@ #define CABANA_VERLETLIST_HPP #include +#include #include #include #include @@ -85,24 +86,40 @@ struct LinkedCellStencil template struct VerletListBuilder + : public N2NeighborListBuilder { // Types. + using n2_list_type = + N2NeighborListBuilder; + using PositionValueType = typename n2_list_type::PositionValueType; using device = DeviceType; - using PositionValueType = typename PositionSlice::value_type; - using RandomAccessPositionSlice = - typename PositionSlice::random_access_slice; using memory_space = typename device::memory_space; using execution_space = typename device::execution_space; // List data. - NeighborListData _data; + using n2_list_type::_data; // Neighbor cutoff. - PositionValueType rsqr; + using n2_list_type::rsqr; // Positions. - RandomAccessPositionSlice position; - std::size_t pid_begin, pid_end; + using n2_list_type::pid_begin; + using n2_list_type::pid_end; + using n2_list_type::position; + + // Check to count or refill. + using n2_list_type::count; + using n2_list_type::refill; + + // Maximum neighbors per particle + using n2_list_type::max_n; + + using CountNeighborsPolicy = typename n2_list_type::CountNeighborsPolicy; + using CountNeighborsTag = typename n2_list_type::CountNeighborsTag; + using FillNeighborsPolicy = typename n2_list_type::FillNeighborsPolicy; + using FillNeighborsTag = typename n2_list_type::FillNeighborsTag; // Binning Data. BinningData bin_data_1d; @@ -111,13 +128,6 @@ struct VerletListBuilder // Cell stencil. LinkedCellStencil cell_stencil; - // Check to count or refill. - bool refill; - bool count; - - // Maximum neighbors per particle - std::size_t max_n; - // Constructor. VerletListBuilder( PositionSlice slice, const std::size_t begin, const std::size_t end, @@ -126,25 +136,10 @@ struct VerletListBuilder const PositionValueType grid_min[3], const PositionValueType grid_max[3], const std::size_t max_neigh ) - : pid_begin( begin ) - , pid_end( end ) + : n2_list_type( slice, begin, end, neighborhood_radius, max_neigh ) , cell_stencil( neighborhood_radius, cell_size_ratio, grid_min, grid_max ) - , max_n( max_neigh ) { - count = true; - refill = false; - - // Create the count view. - _data.counts = - Kokkos::View( "num_neighbors", slice.size() ); - - // Make a guess for the number of neighbors per particle for 2D lists. - initCounts( LayoutTag() ); - - // Get the positions with random access read-only memory. - position = slice; - // Bin the particles in the grid. Don't actually sort them but make a // permutation vector. Note that we are binning all particles here and // not just the requested range. This is because all particles are @@ -154,20 +149,9 @@ struct VerletListBuilder linked_cell_list = LinkedCellList( position, grid_delta, grid_min, grid_max ); bin_data_1d = linked_cell_list.binningData(); - - // We will use the square of the distance for neighbor determination. - rsqr = neighborhood_radius * neighborhood_radius; } // Neighbor count team operator (only used for CSR lists). - struct CountNeighborsTag - { - }; - using CountNeighborsPolicy = - Kokkos::TeamPolicy, - Kokkos::Schedule>; - KOKKOS_INLINE_FUNCTION void operator()( const CountNeighborsTag&, @@ -239,8 +223,12 @@ struct VerletListBuilder { Kokkos::parallel_reduce( Kokkos::ThreadVectorRange( team, num_n ), - [&]( const int n, int& local_count ) { - neighbor_kernel( pid, x_p, y_p, z_p, n_offset, n, local_count ); + [&]( const int n, int& local_count ) + { + // Get the true id of the candidate neighbor. + std::size_t nid = linked_cell_list.permutation( n_offset + n ); + + this->neighbor_kernel( pid, x_p, y_p, z_p, nid, local_count ); }, cell_count ); } @@ -254,136 +242,14 @@ struct VerletListBuilder TeamOpTag ) const { for ( int n = 0; n < num_n; n++ ) - neighbor_kernel( pid, x_p, y_p, z_p, n_offset, n, cell_count ); - } - - // Neighbor count kernel - KOKKOS_INLINE_FUNCTION - void neighbor_kernel( const int pid, const double x_p, const double y_p, - const double z_p, const int n_offset, const int n, - int& local_count ) const - { - // Get the true id of the candidate neighbor. - std::size_t nid = linked_cell_list.permutation( n_offset + n ); - - // Cache the candidate neighbor particle coordinates. - double x_n = position( nid, 0 ); - double y_n = position( nid, 1 ); - double z_n = position( nid, 2 ); - - // If this could be a valid neighbor, continue. - if ( NeighborDiscriminator::isValid( - pid, x_p, y_p, z_p, nid, x_n, y_n, z_n ) ) - { - // Calculate the distance between the particle and its candidate - // neighbor. - PositionValueType dx = x_p - x_n; - PositionValueType dy = y_p - y_n; - PositionValueType dz = z_p - z_n; - PositionValueType dist_sqr = dx * dx + dy * dy + dz * dz; - - // If within the cutoff add to the count. - if ( dist_sqr <= rsqr ) - local_count += 1; - } - } - - // Process the CSR counts by computing offsets and allocating the neighbor - // list. - template - struct OffsetScanOp - { - using kokkos_mem_space = KokkosMemorySpace; - Kokkos::View counts; - Kokkos::View offsets; - KOKKOS_INLINE_FUNCTION - void operator()( const int i, int& update, const bool final_pass ) const - { - if ( final_pass ) - offsets( i ) = update; - update += counts( i ); - } - }; - - void initCounts( VerletLayoutCSR ) {} - - void initCounts( VerletLayout2D ) - { - if ( max_n > 0 ) { - count = false; - - _data.neighbors = Kokkos::View( - Kokkos::ViewAllocateWithoutInitializing( "neighbors" ), - _data.counts.size(), max_n ); - } - } - - void processCounts( VerletLayoutCSR ) - { - // Allocate offsets. - _data.offsets = Kokkos::View( - Kokkos::ViewAllocateWithoutInitializing( "neighbor_offsets" ), - _data.counts.size() ); - - // Calculate offsets from counts and the total number of counts. - OffsetScanOp offset_op; - offset_op.counts = _data.counts; - offset_op.offsets = _data.offsets; - int total_num_neighbor; - Kokkos::RangePolicy range_policy( - 0, _data.counts.extent( 0 ) ); - Kokkos::parallel_scan( "Cabana::VerletListBuilder::offset_scan", - range_policy, offset_op, total_num_neighbor ); - Kokkos::fence(); - - // Allocate the neighbor list. - _data.neighbors = Kokkos::View( - Kokkos::ViewAllocateWithoutInitializing( "neighbors" ), - total_num_neighbor ); - - // Reset the counts. We count again when we fill. - Kokkos::deep_copy( _data.counts, 0 ); - } - - // Process 2D counts by computing the maximum number of neighbors and - // reallocating the 2D data structure if needed. - void processCounts( VerletLayout2D ) - { - // Calculate the maximum number of neighbors. - auto counts = _data.counts; - int max_num_neighbor = 0; - Kokkos::Max max_reduce( max_num_neighbor ); - Kokkos::parallel_reduce( - "Cabana::VerletListBuilder::reduce_max", - Kokkos::RangePolicy( 0, _data.counts.size() ), - KOKKOS_LAMBDA( const int i, int& value ) { - if ( counts( i ) > value ) - value = counts( i ); - }, - max_reduce ); - Kokkos::fence(); + // Get the true id of the candidate neighbor. + std::size_t nid = linked_cell_list.permutation( n_offset + n ); - // Reallocate the neighbor list if previous size is exceeded. - if ( count or ( std::size_t ) - max_num_neighbor > _data.neighbors.extent( 1 ) ) - { - refill = true; - Kokkos::deep_copy( _data.counts, 0 ); - _data.neighbors = Kokkos::View( - Kokkos::ViewAllocateWithoutInitializing( "neighbors" ), - _data.counts.size(), max_num_neighbor ); + this->neighbor_kernel( pid, x_p, y_p, z_p, nid, cell_count ); } } - // Neighbor count team operator. - struct FillNeighborsTag - { - }; - using FillNeighborsPolicy = - Kokkos::TeamPolicy, - Kokkos::Schedule>; KOKKOS_INLINE_FUNCTION void operator()( const FillNeighborsTag&, @@ -447,8 +313,14 @@ struct VerletListBuilder TeamVectorOpTag ) const { Kokkos::parallel_for( - Kokkos::ThreadVectorRange( team, num_n ), [&]( const int n ) - { neighbor_kernel( pid, x_p, y_p, z_p, n_offset, n ); } ); + Kokkos::ThreadVectorRange( team, num_n ), + [&]( const int n ) + { + // Get the true id of the candidate neighbor. + std::size_t nid = linked_cell_list.permutation( n_offset + n ); + + this->neighbor_kernel( pid, x_p, y_p, z_p, nid ); + } ); } // Neighbor fill serial loop. @@ -461,41 +333,14 @@ struct VerletListBuilder for ( int n = 0; n < num_n; n++ ) Kokkos::single( Kokkos::PerThread( team ), - [&]() { neighbor_kernel( pid, x_p, y_p, z_p, n_offset, n ); } ); - } - - // Neighbor fill kernel. - KOKKOS_INLINE_FUNCTION - void neighbor_kernel( const int pid, const double x_p, const double y_p, - const double z_p, const int n_offset, - const int n ) const - { - // Get the true id of the candidate neighbor. - std::size_t nid = linked_cell_list.permutation( n_offset + n ); - - // Cache the candidate neighbor particle coordinates. - double x_n = position( nid, 0 ); - double y_n = position( nid, 1 ); - double z_n = position( nid, 2 ); + [&]() + { + // Get the true id of the candidate neighbor. + std::size_t nid = + linked_cell_list.permutation( n_offset + n ); - // If this could be a valid neighbor, continue. - if ( NeighborDiscriminator::isValid( - pid, x_p, y_p, z_p, nid, x_n, y_n, z_n ) ) - { - // Calculate the distance between the particle and its candidate - // neighbor. - PositionValueType dx = x_p - x_n; - PositionValueType dy = y_p - y_n; - PositionValueType dz = z_p - z_n; - PositionValueType dist_sqr = dx * dx + dy * dy + dz * dz; - - // If within the cutoff increment the neighbor count and add as a - // neighbor at that index. - if ( dist_sqr <= rsqr ) - { - _data.addNeighbor( pid, nid ); - } - } + this->neighbor_kernel( pid, x_p, y_p, z_p, nid ); + } ); } }; diff --git a/core/unit_test/CMakeLists.txt b/core/unit_test/CMakeLists.txt index eac64bfcf..969855c94 100644 --- a/core/unit_test/CMakeLists.txt +++ b/core/unit_test/CMakeLists.txt @@ -22,6 +22,7 @@ set(SERIAL_TESTS AoSoA DeepCopy LinkedCellList + NeighborListN2 NeighborList Parallel ParameterPack diff --git a/core/unit_test/tstNeighborList.hpp b/core/unit_test/tstNeighborList.hpp index d00fde467..c542e1b0b 100644 --- a/core/unit_test/tstNeighborList.hpp +++ b/core/unit_test/tstNeighborList.hpp @@ -255,16 +255,17 @@ TEST( TEST_CATEGORY, verlet_list_half_test ) TEST( TEST_CATEGORY, verlet_list_full_range_test ) { #ifndef KOKKOS_ENABLE_OPENMPTARGET // FIXME_OPENMPTARGET - testVerletListFullPartialRange(); #endif - testVerletListFullPartialRange(); + testVerletListFullPartialRange(); #ifndef KOKKOS_ENABLE_OPENMPTARGET // FIXME_OPENMPTARGET - testVerletListFullPartialRange(); #endif - testVerletListFullPartialRange(); } diff --git a/core/unit_test/tstNeighborListN2.hpp b/core/unit_test/tstNeighborListN2.hpp new file mode 100644 index 000000000..610ce2172 --- /dev/null +++ b/core/unit_test/tstNeighborListN2.hpp @@ -0,0 +1,201 @@ +/**************************************************************************** + * Copyright (c) 2018-2022 by the Cabana authors * + * All rights reserved. * + * * + * This file is part of the Cabana library. Cabana is distributed under a * + * BSD 3-clause license. For the licensing terms see the LICENSE file in * + * the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +#include +#include +#include +#include + +#include + +#include + +#include + +namespace Test +{ + +//---------------------------------------------------------------------------// +template +void testN2List() +{ + // Create the AoSoA and fill with random particle positions. + NeighborListTestData test_data; + auto position = Cabana::slice<0>( test_data.aosoa ); + + // Create the neighbor list. + { + Cabana::N2NeighborList + nlist_full( position, 0, position.size(), test_data.test_radius ); + // Test default construction. + Cabana::N2NeighborList + nlist; + + nlist = nlist_full; + + checkNeighborList( nlist, test_data.N2_list_copy, + test_data.num_particle, Cabana::FullNeighborTag{} ); + + // Test rebuild function with explict execution space. + nlist.build( TEST_EXECSPACE{}, position, 0, position.size(), + test_data.test_radius ); + checkNeighborList( nlist, test_data.N2_list_copy, + test_data.num_particle, Cabana::FullNeighborTag{} ); + } + // Check again, building with a large array allocation size + { + Cabana::N2NeighborList + nlist_max( position, 0, position.size(), test_data.test_radius, + 100 ); + checkNeighborList( nlist_max, test_data.N2_list_copy, + test_data.num_particle, Cabana::FullNeighborTag{} ); + } + // Check again, building with a small array allocation size (refill) + { + Cabana::N2NeighborList + nlist_max2( position, 0, position.size(), test_data.test_radius, + 2 ); + checkNeighborList( nlist_max2, test_data.N2_list_copy, + test_data.num_particle, Cabana::FullNeighborTag{} ); + } +} + +//---------------------------------------------------------------------------// +template +void testN2ListFullPartialRange() +{ + // Create the AoSoA and fill with random particle positions. + NeighborListTestData test_data; + auto position = Cabana::slice<0>( test_data.aosoa ); + + // Create the neighbor list. + Cabana::N2NeighborList + nlist( position, 0, test_data.num_ignore, test_data.test_radius ); + + // Check the neighbor list. + checkFullNeighborListPartialRange( nlist, test_data.N2_list_copy, + test_data.num_particle, + test_data.num_ignore ); +} + +template +void testN2NeighborParallelFor() +{ + // Create the AoSoA and fill with random particle positions. + NeighborListTestData test_data; + auto position = Cabana::slice<0>( test_data.aosoa ); + + // Create the neighbor list. + using ListType = + Cabana::N2NeighborList; + ListType nlist( position, 0, position.size(), test_data.test_radius ); + + checkNeighborParallelFor( nlist, test_data.N2_list_copy, + test_data.num_particle ); +} + +template +void testN2NeighborParallelReduce() +{ + // Create the AoSoA and fill with random particle positions. + NeighborListTestData test_data; + auto position = Cabana::slice<0>( test_data.aosoa ); + + // Create the neighbor list. + using ListType = + Cabana::N2NeighborList; + ListType nlist( position, 0, position.size(), test_data.test_radius ); + + checkNeighborParallelReduce( nlist, test_data.N2_list_copy, + test_data.aosoa ); +} + +//---------------------------------------------------------------------------// +// TESTS +//---------------------------------------------------------------------------// +TEST( TEST_CATEGORY, n2_list_full_test ) +{ +#ifndef KOKKOS_ENABLE_OPENMPTARGET // FIXME_OPENMPTARGET + testN2List(); +#endif + testN2List(); + +#ifndef KOKKOS_ENABLE_OPENMPTARGET // FIXME_OPENMPTARGET + testN2List(); +#endif + testN2List(); +} + +//---------------------------------------------------------------------------// +TEST( TEST_CATEGORY, n2_list_half_test ) +{ +#ifndef KOKKOS_ENABLE_OPENMPTARGET // FIXME_OPENMPTARGET + testN2List(); +#endif + testN2List(); + +#ifndef KOKKOS_ENABLE_OPENMPTARGET // FIXME_OPENMPTARGET + testN2List(); +#endif + testN2List(); +} + +//---------------------------------------------------------------------------// +TEST( TEST_CATEGORY, n2_list_full_range_test ) +{ +#ifndef KOKKOS_ENABLE_OPENMPTARGET // FIXME_OPENMPTARGET + testN2ListFullPartialRange(); +#endif + testN2ListFullPartialRange(); + +#ifndef KOKKOS_ENABLE_OPENMPTARGET // FIXME_OPENMPTARGET + testN2ListFullPartialRange(); +#endif + testN2ListFullPartialRange(); +} + +//---------------------------------------------------------------------------// +TEST( TEST_CATEGORY, n2_parallel_for_test ) +{ +#ifndef KOKKOS_ENABLE_OPENMPTARGET // FIXME_OPENMPTARGET + testN2NeighborParallelFor(); +#endif + testN2NeighborParallelFor(); +} + +//---------------------------------------------------------------------------// +TEST( TEST_CATEGORY, n2_parallel_reduce_test ) +{ +#ifndef KOKKOS_ENABLE_OPENMPTARGET // FIXME_OPENMPTARGET + testN2NeighborParallelReduce(); +#endif + testN2NeighborParallelReduce(); +} +//---------------------------------------------------------------------------// + +} // namespace Test