diff --git a/cajita/src/CMakeLists.txt b/cajita/src/CMakeLists.txt index 2315321b5..16d0f532d 100644 --- a/cajita/src/CMakeLists.txt +++ b/cajita/src/CMakeLists.txt @@ -29,6 +29,7 @@ set(HEADERS_PUBLIC Cajita_ManualPartitioner.hpp Cajita_MpiTraits.hpp Cajita_Parallel.hpp + Cajita_ParticleGridDistributor.hpp Cajita_Partitioner.hpp Cajita_ReferenceStructuredSolver.hpp Cajita_Splines.hpp diff --git a/cajita/src/Cajita.hpp b/cajita/src/Cajita.hpp index 5525df5f8..8ccee93ef 100644 --- a/cajita/src/Cajita.hpp +++ b/cajita/src/Cajita.hpp @@ -31,6 +31,7 @@ #include #include #include +#include #include #include #include diff --git a/cajita/src/Cajita_ParticleGridDistributor.hpp b/cajita/src/Cajita_ParticleGridDistributor.hpp new file mode 100644 index 000000000..542abc4d1 --- /dev/null +++ b/cajita/src/Cajita_ParticleGridDistributor.hpp @@ -0,0 +1,375 @@ +/**************************************************************************** + * Copyright (c) 2018-2021 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 * + ****************************************************************************/ + +#ifndef CABANA_PARTICLEGRIDDISTRIBUTOR_HPP +#define CABANA_PARTICLEGRIDDISTRIBUTOR_HPP + +#include +#include + +#include +#include +#include +#include + +#include + +#include + +namespace Cajita +{ + +//---------------------------------------------------------------------------// +// Particle Grid Distributor +//---------------------------------------------------------------------------// + +namespace Impl +{ +//! \cond Impl + +// Build neighbor topology of 27 nearest 3D neighbors. Some of the ranks in +// this list may be invalid. +template +std::enable_if_t<3 == NSD, std::vector> +getTopology( const LocalGridType& local_grid ) +{ + std::vector topology( 27, -1 ); + int nr = 0; + for ( int k = -1; k < 2; ++k ) + for ( int j = -1; j < 2; ++j ) + for ( int i = -1; i < 2; ++i, ++nr ) + topology[nr] = local_grid.neighborRank( i, j, k ); + return topology; +} + +// Build neighbor topology of 8 nearest 2D neighbors. Some of the ranks in +// this list may be invalid. +template +std::enable_if_t<2 == NSD, std::vector> +getTopology( const LocalGridType& local_grid ) +{ + std::vector topology( 9, -1 ); + int nr = 0; + for ( int j = -1; j < 2; ++j ) + for ( int i = -1; i < 2; ++i, ++nr ) + topology[nr] = local_grid.neighborRank( i, j ); + return topology; +} + +// Locate the particles in the local grid and get their destination rank. +// Particles are assumed to only migrate to a location in the nearest +// neighbor halo or stay on this rank. If the particle crosses a global +// periodic boundary, wrap it's coordinates back into the domain. +template +void getMigrateDestinations( const LocalGridType& local_grid, + const NeighborRankView& neighbor_ranks, + DestinationRankView& destinations, + PositionSliceType& positions ) +{ + static constexpr std::size_t num_space_dim = LocalGridType::num_space_dim; + using execution_space = typename PositionSliceType::execution_space; + + // Check within the local domain. + const auto& local_mesh = + Cajita::createLocalMesh( local_grid ); + + // Use global domain for periodicity. + const auto& global_grid = local_grid.globalGrid(); + const auto& global_mesh = global_grid.globalMesh(); + + Kokkos::Array local_low{}; + Kokkos::Array local_high{}; + Kokkos::Array periodic{}; + Kokkos::Array global_low{}; + Kokkos::Array global_high{}; + Kokkos::Array global_extent{}; + + for ( std::size_t d = 0; d < num_space_dim; ++d ) + { + local_low[d] = local_mesh.lowCorner( Cajita::Own(), d ); + local_high[d] = local_mesh.highCorner( Cajita::Own(), d ); + periodic[d] = global_grid.isPeriodic( d ); + global_low[d] = global_mesh.lowCorner( d ); + global_high[d] = global_mesh.highCorner( d ); + global_extent[d] = global_mesh.extent( d ); + } + + Kokkos::parallel_for( + "get_migrate_destinations", + Kokkos::RangePolicy( 0, positions.size() ), + KOKKOS_LAMBDA( const int p ) { + // Compute the logical index of the neighbor we are sending to. + int nid[num_space_dim]; + for ( std::size_t d = 0; d < num_space_dim; ++d ) + { + nid[d] = 1; + if ( positions( p, d ) < local_low[d] ) + nid[d] = 0; + else if ( positions( p, d ) > local_high[d] ) + nid[d] = 2; + } + + // Compute the destination MPI rank [ni + 3*(nj + 3*nk) in 3d]. + int neighbor_index = nid[0]; + for ( std::size_t d = 1; d < num_space_dim; ++d ) + { + int npower = 1; + for ( std::size_t dp = 1; dp <= d; ++dp ) + npower *= 3; + neighbor_index += npower * nid[d]; + } + destinations( p ) = neighbor_ranks( neighbor_index ); + + // Shift particles through periodic boundaries. + for ( std::size_t d = 0; d < num_space_dim; ++d ) + { + if ( periodic[d] ) + { + if ( positions( p, d ) > global_high[d] ) + positions( p, d ) -= global_extent[d]; + else if ( positions( p, d ) < global_low[d] ) + positions( p, d ) += global_extent[d]; + } + } + } ); +} +//! \endcond +} // namespace Impl + +//-----------------------------------------------------------------------// +/*! + \brief Check for the number of particles that must be communicated + + \tparam LocalGridType Cajita LocalGrid type. + + \tparam PositionSliceType Particle position type. + + \param local_grid The local grid containing periodicity and system bound + information. + + \param positions The particle position container, either Slice or View. + + \param minimum_halo_width Number of halo mesh widths to include for + ghosting. +*/ +template +int migrateCount( const LocalGridType& local_grid, + const PositionSliceType& positions, + const int minimum_halo_width ) +{ + using grid_type = LocalGridType; + static constexpr std::size_t num_space_dim = grid_type::num_space_dim; + using mesh_type = typename grid_type::mesh_type; + using scalar_type = typename mesh_type::scalar_type; + using uniform_type = UniformMesh; + static_assert( std::is_same::value, + "Migrate count requires a uniform mesh." ); + + using execution_space = typename PositionSliceType::execution_space; + + // Check within the halo width, within the ghosted domain. + const auto& local_mesh = + Cajita::createLocalMesh( local_grid ); + + Kokkos::Array local_low{}; + Kokkos::Array local_high{}; + for ( std::size_t d = 0; d < num_space_dim; ++d ) + { + auto dx = local_grid.globalGrid().globalMesh().cellSize( d ); + local_low[d] = local_mesh.lowCorner( Cajita::Ghost(), d ) + + minimum_halo_width * dx; + local_high[d] = local_mesh.highCorner( Cajita::Ghost(), d ) - + minimum_halo_width * dx; + } + + int comm_count = 0; + Kokkos::parallel_reduce( + "redistribute_count", + Kokkos::RangePolicy( 0, positions.size() ), + KOKKOS_LAMBDA( const int p, int& result ) { + for ( std::size_t d = 0; d < num_space_dim; ++d ) + if ( positions( p, d ) < local_low[d] || + positions( p, d ) > local_high[d] ) + { + result += 1; + break; + } + }, + comm_count ); + + MPI_Allreduce( MPI_IN_PLACE, &comm_count, 1, MPI_INT, MPI_SUM, + local_grid.globalGrid().comm() ); + + return comm_count; +} + +//---------------------------------------------------------------------------// +/*! + \brief Determine which data should be migrated from one uniquely-owned + decomposition to another uniquely-owned decomposition, using bounds of a + Cajita grid and taking periodic boundaries into account. + + \tparam LocalGridType Cajita LocalGrid type. + + \tparam PositionSliceType Position type. + + \param local_grid The local grid containing periodicity and system bound + information. + + \param positions The particle positions. + + \return Distributor for later migration. +*/ +template +Cabana::Distributor +createParticleGridDistributor( const LocalGridType& local_grid, + PositionSliceType& positions ) +{ + using device_type = typename PositionSliceType::device_type; + + // Get all 26 neighbor ranks. + auto topology = Impl::getTopology( local_grid ); + + Kokkos::View + topology_host( topology.data(), topology.size() ); + auto topology_mirror = + Kokkos::create_mirror_view_and_copy( device_type(), topology_host ); + Kokkos::View destinations( + Kokkos::ViewAllocateWithoutInitializing( "destinations" ), + positions.size() ); + + // Determine destination ranks for all particles and wrap positions across + // periodic boundaries. + Impl::getMigrateDestinations( local_grid, topology_mirror, destinations, + positions ); + + // Create the Cabana distributor. + Cabana::Distributor distributor( + local_grid.globalGrid().comm(), destinations, topology ); + return distributor; +} + +//---------------------------------------------------------------------------// +/*! + \brief Migrate data from one uniquely-owned decomposition to another + uniquely-owned decomposition, using the bounds and periodic boundaries of a + Cajita grid to determine which particles should be moved. In-place variant. + + \tparam LocalGridType Cajita LocalGrid type. + + \tparam ParticlePositions Particle position type. + + \tparam PositionContainer AoSoA type. + + \param local_grid The local grid containing periodicity and system bounds. + + \param positions Particle positions. + + \param particles The particle AoSoA. + + \param min_halo_width Number of halo mesh widths to allow particles before + migrating. + + \param force_migrate Migrate particles outside the local domain regardless of + ghosted halo. +*/ +template +void particleGridMigrate( const LocalGridType& local_grid, + ParticlePositions& positions, + ParticleContainer& particles, + const int min_halo_width, + const bool force_migrate = false ) +{ + // When false, this option checks that any particles are nearly outside the + // ghosted halo region (outside the min_halo_width) before initiating + // migration. Otherwise, anything outside the local domain is migrated + // regardless of position in the halo. + if ( !force_migrate ) + { + // Check to see if we need to communicate. + auto comm_count = migrateCount( local_grid, positions, min_halo_width ); + + // If we have no particles near the ghosted boundary, then exit. + if ( 0 == comm_count ) + return; + } + + auto distributor = createParticleGridDistributor( local_grid, positions ); + + // Redistribute the particles. + migrate( distributor, particles ); +} + +//---------------------------------------------------------------------------// +/*! + \brief Migrate data from one uniquely-owned decomposition to another + uniquely-owned decomposition, using the bounds and periodic boundaries of a + Cajita grid to determine which particles should be moved. Separate AoSoA + variant. + + \tparam LocalGridType Cajita LocalGrid type. + + \tparam ParticlePositions Particle position type. + + \tparam ParticleContainer AoSoA type. + + \param local_grid The local grid containing periodicity and system bounds. + + \param positions Particle positions. + + \param src_particles The source particle AoSoA. + + \param dst_particles The destination particle AoSoA. + + \param min_halo_width Number of halo mesh widths to allow particles before + migrating. + + \param force_migrate Migrate particles outside the local domain regardless of + ghosted halo. +*/ +template +void particleGridMigrate( const LocalGridType& local_grid, + const ParticlePositions& positions, + const ParticleContainer& src_particles, + ParticleContainer& dst_particles, + const int min_halo_width, + const bool force_migrate = false ) +{ + // When false, this option checks that any particles are nearly outside the + // ghosted halo region (outside the min_halo_width) before initiating + // migration. Otherwise, anything outside the local domain is migrated + // regardless of position in the halo. + if ( !force_migrate ) + { + // Check to see if we need to communicate. + auto comm_count = migrateCount( local_grid, positions, min_halo_width ); + + // If we have no particles near the ghosted boundary, copy, then exit. + if ( 0 == comm_count ) + { + Cabana::deep_copy( dst_particles, src_particles ); + return; + } + } + + auto distributor = createParticleGridDistributor( local_grid, positions ); + + // Resize as needed. + dst_particles.resize( distributor.totalNumImport() ); + + // Redistribute the particles. + migrate( distributor, src_particles, dst_particles ); +} + +} // namespace Cajita + +#endif // end CABANA_PARTICLEGRIDDISTRIBUTOR_HPP diff --git a/cajita/unit_test/CMakeLists.txt b/cajita/unit_test/CMakeLists.txt index b810c41df..0c831c1ff 100644 --- a/cajita/unit_test/CMakeLists.txt +++ b/cajita/unit_test/CMakeLists.txt @@ -28,6 +28,8 @@ set(MPI_TESTS Array2d Halo3d Halo2d + ParticleGridDistributor2d + ParticleGridDistributor3d SplineEvaluation3d SplineEvaluation2d Interpolation3d diff --git a/cajita/unit_test/tstParticleGridDistributor2d.hpp b/cajita/unit_test/tstParticleGridDistributor2d.hpp new file mode 100644 index 000000000..4d2d7490f --- /dev/null +++ b/cajita/unit_test/tstParticleGridDistributor2d.hpp @@ -0,0 +1,333 @@ +/**************************************************************************** + * Copyright (c) 2018-2021 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 +#include +#include +#include + +#include + +#include + +#include + +#include +#include +#include + +namespace Test +{ + +using Cajita::Dim; + +//---------------------------------------------------------------------------// +template +void migrateTest( const GridType global_grid, const double cell_size, + const int data_halo_size, const int test_halo_size, + const bool force_comm, const int test_type ) +{ + // Create local block with varying halo size. + auto block = Cajita::createLocalGrid( global_grid, data_halo_size ); + auto local_mesh = Cajita::createLocalMesh( *block ); + + // Allocate a maximum number of particles assuming we have a halo on every + // boundary. + auto ghosted_cell_space = + block->indexSpace( Cajita::Ghost(), Cajita::Cell(), Cajita::Local() ); + int num_particle = ghosted_cell_space.size(); + using MemberTypes = Cabana::MemberTypes; + using ParticleContainer = Cabana::AoSoA; + ParticleContainer particles( "particles", num_particle ); + auto coords = Cabana::slice<0>( particles, "coords" ); + auto linear_ids = Cabana::slice<1>( particles, "linear_ids" ); + + // Put particles in the center of every cell including halo cells if we + // have them. Their ids should be equivalent to that of the rank they are + // going to. + int pid = 0; + for ( int nj = -1; nj < 2; ++nj ) + for ( int ni = -1; ni < 2; ++ni ) + { + auto neighbor_rank = block->neighborRank( ni, nj ); + if ( neighbor_rank >= 0 ) + { + auto shared_space = block->sharedIndexSpace( + Cajita::Ghost(), Cajita::Cell(), ni, nj ); + for ( int j = shared_space.min( Dim::J ); + j < shared_space.max( Dim::J ); ++j ) + for ( int i = shared_space.min( Dim::I ); + i < shared_space.max( Dim::I ); ++i ) + { + // Set the coordinates at the cell center. + coords( pid, Dim::I ) = + local_mesh.lowCorner( Cajita::Ghost(), Dim::I ) + + ( i + 0.5 ) * cell_size; + coords( pid, Dim::J ) = + local_mesh.lowCorner( Cajita::Ghost(), Dim::J ) + + ( j + 0.5 ) * cell_size; + + // Set the linear ids as the linear rank of + // the neighbor. + linear_ids( pid ) = neighbor_rank; + + // Increment the particle count. + ++pid; + } + } + } + num_particle = pid; + particles.resize( num_particle ); + + ParticleContainer particles_initial( "initial_particles", num_particle ); + Cabana::deep_copy( particles_initial, particles ); + auto coords_initial = Cabana::slice<0>( particles_initial ); + + // Copy to the device space. + auto particles_mirror = + Cabana::create_mirror_view_and_copy( TEST_DEVICE(), particles ); + auto coords_mirror = Cabana::slice<0>( particles_mirror, "coords" ); + + // Redistribute the particle AoSoA in place. + if ( test_type == 0 ) + { + Cajita::particleGridMigrate( *block, coords_mirror, particles_mirror, + test_halo_size, force_comm ); + + // Copy back to check. + particles = Cabana::create_mirror_view_and_copy( Kokkos::HostSpace(), + particles_mirror ); + } + // Do the migration with a separate destination AoSoA. + else if ( test_type == 1 ) + { + auto particles_dst = + Cabana::create_mirror_view( TEST_MEMSPACE(), particles_mirror ); + Cajita::particleGridMigrate( *block, coords_mirror, particles_mirror, + particles_dst, test_halo_size, + force_comm ); + // Copy back to check. + particles = Cabana::create_mirror_view_and_copy( Kokkos::HostSpace(), + particles_dst ); + } + + coords = Cabana::slice<0>( particles, "coords" ); + linear_ids = Cabana::slice<1>( particles, "linear_ids" ); + + // Check that we got as many particles as we should have. + EXPECT_EQ( coords.size(), num_particle ); + EXPECT_EQ( linear_ids.size(), num_particle ); + + std::array local_low = { + local_mesh.lowCorner( Cajita::Own(), Dim::I ), + local_mesh.lowCorner( Cajita::Own(), Dim::J ) }; + std::array local_high = { + local_mesh.highCorner( Cajita::Own(), Dim::I ), + local_mesh.highCorner( Cajita::Own(), Dim::J ) }; + + for ( int p = 0; p < num_particle; ++p ) + { + // Particles should be redistributed if forcing communication or if + // anything is outside the minimum halo width (currently every case + // except test_halo_size=0) + if ( force_comm || test_halo_size > 0 ) + { + // Check that all of the particle ids are equal to this rank id. + EXPECT_EQ( linear_ids( p ), global_grid->blockId() ); + + // Check that all of the particles are now in the local domain. + for ( int d = 0; d < 2; ++d ) + { + EXPECT_GE( coords( p, d ), local_low[d] ); + EXPECT_LE( coords( p, d ), local_high[d] ); + } + } + else + { + // If nothing was outside the halo, nothing should have moved. + for ( int d = 0; d < 2; ++d ) + EXPECT_DOUBLE_EQ( coords( p, d ), coords_initial( p, d ) ); + } + } +} // namespace Test + +//---------------------------------------------------------------------------// +// The objective of this test is to check how the redistribution works when we +// have no particles to redistribute. In this case we put no particles in the +// halo so no communication should occur. This ensures the graph communication +// works when some neighbors get no data. +template +void localOnlyTest( const GridType global_grid, const double cell_size ) +{ + // Get the local block with a halo of 2. + const int data_halo_size = 2; + auto block = Cajita::createLocalGrid( global_grid, data_halo_size ); + auto local_mesh = Cajita::createLocalMesh( *block ); + + // Allocate particles + auto owned_cell_space = + block->indexSpace( Cajita::Own(), Cajita::Cell(), Cajita::Local() ); + int num_particle = owned_cell_space.size(); + using MemberTypes = Cabana::MemberTypes; + using ParticleContainer = Cabana::AoSoA; + ParticleContainer particles( "particles", num_particle ); + auto coords = Cabana::slice<0>( particles, "coords" ); + auto linear_ids = Cabana::slice<1>( particles, "linear_ids" ); + + // Put particles in the center of every local cell. + int pid = 0; + for ( int j = 0; j < owned_cell_space.extent( Dim::J ); ++j ) + for ( int i = 0; i < owned_cell_space.extent( Dim::I ); ++i ) + { + // Set the coordinates at the cell center. + coords( pid, Dim::I ) = + local_mesh.lowCorner( Cajita::Own(), Dim::I ) + + ( i + 0.5 ) * cell_size; + coords( pid, Dim::J ) = + local_mesh.lowCorner( Cajita::Own(), Dim::J ) + + ( j + 0.5 ) * cell_size; + + // Set the linear rank + linear_ids( pid ) = global_grid->blockId(); + + // Increment the particle count. + ++pid; + } + + // Copy to the device space. + auto particles_mirror = + Cabana::create_mirror_view_and_copy( TEST_DEVICE(), particles ); + + // Redistribute the particles. + auto coords_mirror = Cabana::slice<0>( particles_mirror, "coords" ); + Cajita::particleGridMigrate( *block, coords_mirror, particles_mirror, 0, + true ); + + // Copy back to check. + particles = Cabana::create_mirror_view_and_copy( Kokkos::HostSpace(), + particles_mirror ); + coords = Cabana::slice<0>( particles, "coords" ); + linear_ids = Cabana::slice<1>( particles, "linear_ids" ); + + // Check that we got as many particles as we should have. + EXPECT_EQ( coords.size(), num_particle ); + EXPECT_EQ( linear_ids.size(), num_particle ); + + // Check that all of the particle ids are equal to this rank id. + for ( int p = 0; p < num_particle; ++p ) + EXPECT_EQ( linear_ids( p ), global_grid->blockId() ); + + // Check that all of the particles are now in the local domain. + double local_low[3] = { local_mesh.lowCorner( Cajita::Own(), Dim::I ), + local_mesh.lowCorner( Cajita::Own(), Dim::J ) }; + double local_high[3] = { local_mesh.highCorner( Cajita::Own(), Dim::I ), + local_mesh.highCorner( Cajita::Own(), Dim::J ) }; + for ( int p = 0; p < num_particle; ++p ) + for ( int d = 0; d < 2; ++d ) + { + EXPECT_GE( coords( p, d ), local_low[d] ); + EXPECT_LE( coords( p, d ), local_high[d] ); + } +} + +auto createGrid( const Cajita::ManualBlockPartitioner<2>& partitioner, + const std::array& is_periodic, + const double cell_size ) +{ + // Create the global grid. + std::array global_num_cell = { 24, 17 }; + std::array global_low_corner = { 1.2, 3.3 }; + std::array global_high_corner = { + global_low_corner[0] + cell_size * global_num_cell[0], + global_low_corner[1] + cell_size * global_num_cell[1] }; + auto global_mesh = Cajita::createUniformGlobalMesh( + global_low_corner, global_high_corner, global_num_cell ); + auto global_grid = Cajita::createGlobalGrid( MPI_COMM_WORLD, global_mesh, + is_periodic, partitioner ); + return global_grid; +} + +void testParticleGridMigrate( const bool periodic ) +{ + // Let MPI compute the partitioning for this test. + int comm_size; + MPI_Comm_size( MPI_COMM_WORLD, &comm_size ); + std::array ranks_per_dim = { 0, 0 }; + MPI_Dims_create( comm_size, 2, ranks_per_dim.data() ); + Cajita::ManualBlockPartitioner<2> partitioner( ranks_per_dim ); + + // Boundaries are not periodic. + std::array is_periodic = { periodic, periodic }; + + // Create global grid. + double cell_size = 0.23; + auto global_grid = createGrid( partitioner, is_periodic, cell_size ); + + // Test in-place + // Test multiple system halo sizes + for ( int i = 0; i < 3; i++ ) + // Test multiple minimum_halo_width + for ( int j = 0; j < 3; j++ ) + migrateTest( global_grid, cell_size, 1, j, false, 0 ); + + // Retest with separate destination AoSoA. + migrateTest( global_grid, cell_size, 2, 2, true, 1 ); + + // Test with forced communication. + migrateTest( global_grid, cell_size, 2, 2, true, 0 ); + + // Test with different block configurations to make sure all the + // dimensions get partitioned even at small numbers of ranks. + if ( ranks_per_dim[0] != ranks_per_dim[1] ) + { + std::swap( ranks_per_dim[0], ranks_per_dim[1] ); + partitioner = Cajita::ManualBlockPartitioner<2>( ranks_per_dim ); + migrateTest( global_grid, cell_size, 2, 2, true, 0 ); + } +} + +//---------------------------------------------------------------------------// +// RUN TESTS +//---------------------------------------------------------------------------// +TEST( TEST_CATEGORY, not_periodic_test ) { testParticleGridMigrate( false ); } + +//---------------------------------------------------------------------------// +TEST( TEST_CATEGORY, periodic_test ) { testParticleGridMigrate( true ); } + +//---------------------------------------------------------------------------// +TEST( TEST_CATEGORY, local_only_test ) +{ + // Let MPI compute the partitioning for this test. + int comm_size; + MPI_Comm_size( MPI_COMM_WORLD, &comm_size ); + std::array ranks_per_dim = { 0, 0 }; + MPI_Dims_create( comm_size, 2, ranks_per_dim.data() ); + Cajita::ManualBlockPartitioner<2> partitioner( ranks_per_dim ); + + // Every boundary is periodic + std::array is_periodic = { true, true }; + + // Create global grid. + double cell_size = 0.23; + auto global_grid = createGrid( partitioner, is_periodic, cell_size ); + + localOnlyTest( global_grid, cell_size ); +} +//---------------------------------------------------------------------------// + +} // end namespace Test diff --git a/cajita/unit_test/tstParticleGridDistributor3d.hpp b/cajita/unit_test/tstParticleGridDistributor3d.hpp new file mode 100644 index 000000000..1bde3a268 --- /dev/null +++ b/cajita/unit_test/tstParticleGridDistributor3d.hpp @@ -0,0 +1,363 @@ +/**************************************************************************** + * Copyright (c) 2018-2021 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 +#include +#include +#include + +#include + +#include + +#include + +#include +#include +#include + +namespace Test +{ + +using Cajita::Dim; + +//---------------------------------------------------------------------------// +template +void migrateTest( const GridType global_grid, const double cell_size, + const int data_halo_size, const int test_halo_size, + const bool force_comm, const int test_type ) +{ + // Create local block with varying halo size. + auto block = Cajita::createLocalGrid( global_grid, data_halo_size ); + auto local_mesh = Cajita::createLocalMesh( *block ); + + // Allocate a maximum number of particles assuming we have a halo on every + // boundary. + auto ghosted_cell_space = + block->indexSpace( Cajita::Ghost(), Cajita::Cell(), Cajita::Local() ); + int num_particle = ghosted_cell_space.size(); + using MemberTypes = Cabana::MemberTypes; + using ParticleContainer = Cabana::AoSoA; + ParticleContainer particles( "particles", num_particle ); + auto coords = Cabana::slice<0>( particles, "coords" ); + auto linear_ids = Cabana::slice<1>( particles, "linear_ids" ); + + // Put particles in the center of every cell including halo cells if we + // have them. Their ids should be equivalent to that of the rank they are + // going to. + int pid = 0; + for ( int nk = -1; nk < 2; ++nk ) + for ( int nj = -1; nj < 2; ++nj ) + for ( int ni = -1; ni < 2; ++ni ) + { + auto neighbor_rank = block->neighborRank( ni, nj, nk ); + if ( neighbor_rank >= 0 ) + { + auto shared_space = block->sharedIndexSpace( + Cajita::Ghost(), Cajita::Cell(), ni, nj, nk ); + for ( int k = shared_space.min( Dim::K ); + k < shared_space.max( Dim::K ); ++k ) + for ( int j = shared_space.min( Dim::J ); + j < shared_space.max( Dim::J ); ++j ) + for ( int i = shared_space.min( Dim::I ); + i < shared_space.max( Dim::I ); ++i ) + { + // Set the coordinates at the cell center. + coords( pid, Dim::I ) = + local_mesh.lowCorner( Cajita::Ghost(), + Dim::I ) + + ( i + 0.5 ) * cell_size; + coords( pid, Dim::J ) = + local_mesh.lowCorner( Cajita::Ghost(), + Dim::J ) + + ( j + 0.5 ) * cell_size; + coords( pid, Dim::K ) = + local_mesh.lowCorner( Cajita::Ghost(), + Dim::K ) + + ( k + 0.5 ) * cell_size; + + // Set the linear ids as the linear rank of + // the neighbor. + linear_ids( pid ) = neighbor_rank; + + // Increment the particle count. + ++pid; + } + } + } + num_particle = pid; + particles.resize( num_particle ); + + ParticleContainer particles_initial( "initial_particles", num_particle ); + Cabana::deep_copy( particles_initial, particles ); + auto coords_initial = Cabana::slice<0>( particles_initial ); + + // Copy to the device space. + auto particles_mirror = + Cabana::create_mirror_view_and_copy( TEST_DEVICE(), particles ); + auto coords_mirror = Cabana::slice<0>( particles_mirror, "coords" ); + + // Redistribute the particle AoSoA in place. + if ( test_type == 0 ) + { + Cajita::particleGridMigrate( *block, coords_mirror, particles_mirror, + test_halo_size, force_comm ); + + // Copy back to check. + particles = Cabana::create_mirror_view_and_copy( Kokkos::HostSpace(), + particles_mirror ); + } + // Do the migration with a separate destination AoSoA. + else if ( test_type == 1 ) + { + auto particles_dst = + Cabana::create_mirror_view( TEST_MEMSPACE(), particles_mirror ); + Cajita::particleGridMigrate( *block, coords_mirror, particles_mirror, + particles_dst, test_halo_size, + force_comm ); + // Copy back to check. + particles = Cabana::create_mirror_view_and_copy( Kokkos::HostSpace(), + particles_dst ); + } + + coords = Cabana::slice<0>( particles, "coords" ); + linear_ids = Cabana::slice<1>( particles, "linear_ids" ); + + // Check that we got as many particles as we should have. + EXPECT_EQ( coords.size(), num_particle ); + EXPECT_EQ( linear_ids.size(), num_particle ); + + std::array local_low = { + local_mesh.lowCorner( Cajita::Own(), Dim::I ), + local_mesh.lowCorner( Cajita::Own(), Dim::J ), + local_mesh.lowCorner( Cajita::Own(), Dim::K ) }; + std::array local_high = { + local_mesh.highCorner( Cajita::Own(), Dim::I ), + local_mesh.highCorner( Cajita::Own(), Dim::J ), + local_mesh.highCorner( Cajita::Own(), Dim::K ) }; + + for ( int p = 0; p < num_particle; ++p ) + { + // Particles should be redistributed if forcing communication or if + // anything is outside the minimum halo width (currently every case + // except test_halo_size=0) + if ( force_comm || test_halo_size > 0 ) + { + // Check that all of the particle ids are equal to this rank id. + EXPECT_EQ( linear_ids( p ), global_grid->blockId() ); + + // Check that all of the particles are now in the local domain. + for ( int d = 0; d < 3; ++d ) + { + EXPECT_GE( coords( p, d ), local_low[d] ); + EXPECT_LE( coords( p, d ), local_high[d] ); + } + } + else + { + // If nothing was outside the halo, nothing should have moved. + for ( int d = 0; d < 3; ++d ) + EXPECT_DOUBLE_EQ( coords( p, d ), coords_initial( p, d ) ); + } + } +} // namespace Test + +//---------------------------------------------------------------------------// +// The objective of this test is to check how the redistribution works when we +// have no particles to redistribute. In this case we put no particles in the +// halo so no communication should occur. This ensures the graph communication +// works when some neighbors get no data. +template +void localOnlyTest( const GridType global_grid, const double cell_size ) +{ + // Get the local block with a halo of 2. + const int data_halo_size = 2; + auto block = Cajita::createLocalGrid( global_grid, data_halo_size ); + auto local_mesh = Cajita::createLocalMesh( *block ); + + // Allocate particles + auto owned_cell_space = + block->indexSpace( Cajita::Own(), Cajita::Cell(), Cajita::Local() ); + int num_particle = owned_cell_space.size(); + using MemberTypes = Cabana::MemberTypes; + using ParticleContainer = Cabana::AoSoA; + ParticleContainer particles( "particles", num_particle ); + auto coords = Cabana::slice<0>( particles, "coords" ); + auto linear_ids = Cabana::slice<1>( particles, "linear_ids" ); + + // Put particles in the center of every local cell. + int pid = 0; + for ( int k = 0; k < owned_cell_space.extent( Dim::K ); ++k ) + for ( int j = 0; j < owned_cell_space.extent( Dim::J ); ++j ) + for ( int i = 0; i < owned_cell_space.extent( Dim::I ); ++i ) + { + // Set the coordinates at the cell center. + coords( pid, Dim::I ) = + local_mesh.lowCorner( Cajita::Own(), Dim::I ) + + ( i + 0.5 ) * cell_size; + coords( pid, Dim::J ) = + local_mesh.lowCorner( Cajita::Own(), Dim::J ) + + ( j + 0.5 ) * cell_size; + coords( pid, Dim::K ) = + local_mesh.lowCorner( Cajita::Own(), Dim::K ) + + ( k + 0.5 ) * cell_size; + + // Set the linear rank + linear_ids( pid ) = global_grid->blockId(); + + // Increment the particle count. + ++pid; + } + + // Copy to the device space. + auto particles_mirror = + Cabana::create_mirror_view_and_copy( TEST_DEVICE(), particles ); + + // Redistribute the particles. + auto coords_mirror = Cabana::slice<0>( particles_mirror, "coords" ); + Cajita::particleGridMigrate( *block, coords_mirror, particles_mirror, 0, + true ); + + // Copy back to check. + particles = Cabana::create_mirror_view_and_copy( Kokkos::HostSpace(), + particles_mirror ); + coords = Cabana::slice<0>( particles, "coords" ); + linear_ids = Cabana::slice<1>( particles, "linear_ids" ); + + // Check that we got as many particles as we should have. + EXPECT_EQ( coords.size(), num_particle ); + EXPECT_EQ( linear_ids.size(), num_particle ); + + // Check that all of the particle ids are equal to this rank id. + for ( int p = 0; p < num_particle; ++p ) + EXPECT_EQ( linear_ids( p ), global_grid->blockId() ); + + // Check that all of the particles are now in the local domain. + double local_low[3] = { local_mesh.lowCorner( Cajita::Own(), Dim::I ), + local_mesh.lowCorner( Cajita::Own(), Dim::J ), + local_mesh.lowCorner( Cajita::Own(), Dim::K ) }; + double local_high[3] = { local_mesh.highCorner( Cajita::Own(), Dim::I ), + local_mesh.highCorner( Cajita::Own(), Dim::J ), + local_mesh.highCorner( Cajita::Own(), Dim::K ) }; + for ( int p = 0; p < num_particle; ++p ) + for ( int d = 0; d < 3; ++d ) + { + EXPECT_GE( coords( p, d ), local_low[d] ); + EXPECT_LE( coords( p, d ), local_high[d] ); + } +} + +auto createGrid( const Cajita::ManualPartitioner& partitioner, + const std::array& is_periodic, + const double cell_size ) +{ + // Create the global grid. + std::array global_num_cell = { 18, 15, 9 }; + std::array global_low_corner = { 1.2, 3.3, -2.8 }; + std::array global_high_corner = { + global_low_corner[0] + cell_size * global_num_cell[0], + global_low_corner[1] + cell_size * global_num_cell[1], + global_low_corner[2] + cell_size * global_num_cell[2] }; + auto global_mesh = Cajita::createUniformGlobalMesh( + global_low_corner, global_high_corner, global_num_cell ); + auto global_grid = Cajita::createGlobalGrid( MPI_COMM_WORLD, global_mesh, + is_periodic, partitioner ); + return global_grid; +} + +void testParticleGridMigrate( const bool periodic ) +{ + // Let MPI compute the partitioning for this test. + int comm_size; + MPI_Comm_size( MPI_COMM_WORLD, &comm_size ); + std::array ranks_per_dim = { 0, 0, 0 }; + MPI_Dims_create( comm_size, 3, ranks_per_dim.data() ); + Cajita::ManualPartitioner partitioner( ranks_per_dim ); + + // Every boundary is periodic + std::array is_periodic = { periodic, periodic, periodic }; + + // Create global grid. + double cell_size = 0.23; + auto global_grid = createGrid( partitioner, is_periodic, cell_size ); + + // Test in-place + // Test multiple system halo sizes + for ( int i = 0; i < 3; i++ ) + // Test multiple minimum_halo_width + for ( int j = 0; j < 3; j++ ) + migrateTest( global_grid, cell_size, i, j, false, 0 ); + + // Retest with separate destination AoSoA. + migrateTest( global_grid, cell_size, 2, 2, true, 1 ); + + // Test with forced communication. + migrateTest( global_grid, cell_size, 2, 2, true, 0 ); + + // Test with different block configurations to make sure all the + // dimensions get partitioned even at small numbers of ranks. + if ( ranks_per_dim[0] != ranks_per_dim[1] ) + { + std::swap( ranks_per_dim[0], ranks_per_dim[1] ); + partitioner = Cajita::ManualPartitioner( ranks_per_dim ); + migrateTest( global_grid, cell_size, 2, 2, true, 0 ); + } + if ( ranks_per_dim[0] != ranks_per_dim[2] ) + { + std::swap( ranks_per_dim[0], ranks_per_dim[2] ); + partitioner = Cajita::ManualPartitioner( ranks_per_dim ); + migrateTest( global_grid, cell_size, 2, 2, true, 0 ); + } + if ( ranks_per_dim[1] != ranks_per_dim[2] ) + { + std::swap( ranks_per_dim[1], ranks_per_dim[2] ); + partitioner = Cajita::ManualPartitioner( ranks_per_dim ); + migrateTest( global_grid, cell_size, 2, 2, true, 0 ); + } +} + +//---------------------------------------------------------------------------// +// RUN TESTS +//---------------------------------------------------------------------------// +TEST( TEST_CATEGORY, not_periodic_test ) { testParticleGridMigrate( false ); } + +//---------------------------------------------------------------------------// +TEST( TEST_CATEGORY, periodic_test ) { testParticleGridMigrate( true ); } + +//---------------------------------------------------------------------------// +TEST( TEST_CATEGORY, local_only_test ) +{ + // Let MPI compute the partitioning for this test. + int comm_size; + MPI_Comm_size( MPI_COMM_WORLD, &comm_size ); + std::array ranks_per_dim = { 0, 0, 0 }; + MPI_Dims_create( comm_size, 3, ranks_per_dim.data() ); + Cajita::ManualPartitioner partitioner( ranks_per_dim ); + + // Every boundary is periodic + std::array is_periodic = { true, true, true }; + + // Create global grid. + double cell_size = 0.23; + auto global_grid = createGrid( partitioner, is_periodic, cell_size ); + + localOnlyTest( global_grid, cell_size ); +} +//---------------------------------------------------------------------------// + +} // end namespace Test diff --git a/core/src/Cabana_CommunicationPlan.hpp b/core/src/Cabana_CommunicationPlan.hpp index b6b4ed6f5..ed283e120 100644 --- a/core/src/Cabana_CommunicationPlan.hpp +++ b/core/src/Cabana_CommunicationPlan.hpp @@ -350,6 +350,29 @@ auto countSendsAndCreateSteering( const ExportRankView element_export_ranks, return std::make_pair( neighbor_counts, neighbor_ids ); } +//---------------------------------------------------------------------------// +// Return unique neighbor ranks, with the current rank first. +std::vector getUniqueTopology( std::vector topology ) +{ + auto remove_end = std::remove( topology.begin(), topology.end(), -1 ); + std::sort( topology.begin(), remove_end ); + auto unique_end = std::unique( topology.begin(), remove_end ); + topology.resize( std::distance( topology.begin(), unique_end ) ); + + // Put this rank first. + int my_rank = -1; + MPI_Comm_rank( MPI_COMM_WORLD, &my_rank ); + for ( auto& n : topology ) + { + if ( n == my_rank ) + { + std::swap( n, topology[0] ); + break; + } + } + return topology; +} + //---------------------------------------------------------------------------// //! \endcond } // end namespace Impl @@ -566,11 +589,8 @@ class CommunicationPlan // Store the number of export elements. _num_export_element = element_export_ranks.size(); - // Store the unique neighbors. - _neighbors = neighbor_ranks; - std::sort( _neighbors.begin(), _neighbors.end() ); - auto unique_end = std::unique( _neighbors.begin(), _neighbors.end() ); - _neighbors.resize( std::distance( _neighbors.begin(), unique_end ) ); + // Store the unique neighbors (this rank first). + _neighbors = Impl::getUniqueTopology( neighbor_ranks ); int num_n = _neighbors.size(); // Get the size of this communicator. @@ -585,15 +605,6 @@ class CommunicationPlan // communication space so any mpi tag will do. const int mpi_tag = 1221; - // If we are sending to ourself put that one first in the neighbor - // list. - for ( auto& n : _neighbors ) - if ( n == my_rank ) - { - std::swap( n, _neighbors[0] ); - break; - } - // Initialize import/export sizes. _num_export.assign( num_n, 0 ); _num_import.assign( num_n, 0 ); diff --git a/core/unit_test/tstCommunicationPlan.hpp b/core/unit_test/tstCommunicationPlan.hpp index 5b6c766f7..4396181b4 100644 --- a/core/unit_test/tstCommunicationPlan.hpp +++ b/core/unit_test/tstCommunicationPlan.hpp @@ -576,6 +576,29 @@ void test7( const bool use_topology ) EXPECT_EQ( n, host_steering( n ) ); } +//---------------------------------------------------------------------------// +void testTopology() +{ + // Get my rank. + int my_rank = -1; + MPI_Comm_rank( MPI_COMM_WORLD, &my_rank ); + + std::vector neighbor_ranks = { 0, 2, 2, 1, 1, 1, + 2, 3, 10, 27, my_rank }; + auto unique_ranks = Cabana::Impl::getUniqueTopology( neighbor_ranks ); + + // Check this rank is first. + EXPECT_EQ( my_rank, unique_ranks[0] ); + + // Check each rank is unique. + for ( std::size_t n = 0; n < unique_ranks.size(); ++n ) + for ( std::size_t m = 0; m < unique_ranks.size(); ++m ) + if ( n != m ) + { + EXPECT_NE( unique_ranks[n], unique_ranks[m] ); + } +} + //---------------------------------------------------------------------------// // RUN TESTS //---------------------------------------------------------------------------// @@ -607,6 +630,8 @@ TEST( TEST_CATEGORY, comm_plan_test_6_no_topo ) { test6( false ); } TEST( TEST_CATEGORY, comm_plan_test_7_no_topo ) { test7( false ); } +TEST( TEST_CATEGORY, comm_plan_test_topology ) { testTopology(); } + //---------------------------------------------------------------------------// } // end namespace Test