From 3414152b4566a9ede858cfed3fe3fec358607781 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 20 Feb 2023 14:28:42 -0500 Subject: [PATCH 1/9] Add Picasso grid-based particle initialization Co-authored-by: Stuart Slattery --- cajita/src/Cajita_ParticleInit.hpp | 433 +++++++++++++++++++++++++++ cajita/unit_test/CMakeLists.txt | 1 + cajita/unit_test/tstParticleInit.hpp | 166 ++++++++++ core/src/Cabana_ParticleInit.hpp | 58 ++++ 4 files changed, 658 insertions(+) create mode 100644 cajita/src/Cajita_ParticleInit.hpp create mode 100644 cajita/unit_test/tstParticleInit.hpp diff --git a/cajita/src/Cajita_ParticleInit.hpp b/cajita/src/Cajita_ParticleInit.hpp new file mode 100644 index 000000000..b817535da --- /dev/null +++ b/cajita/src/Cajita_ParticleInit.hpp @@ -0,0 +1,433 @@ +/**************************************************************************** + * 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 * + ****************************************************************************/ + +/**************************************************************************** + * Copyright (c) 2021 by the Picasso authors * + * All rights reserved. * + * * + * This file is part of the Picasso library. Picasso 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 Cajita_ParticleInit.hpp + \brief Particle creation utilities based on uniform grids. +*/ +#ifndef CAJITA_PARTICLEINIT_HPP +#define CAJITA_PARTICLEINIT_HPP + +#include +#include + +#include + +#include +#include + +namespace Cajita +{ + +//---------------------------------------------------------------------------// +/*! + \brief Filter out empty particles that weren't created. + + \param exec_space Kokkos execution space. + \param local_num_create The number of particles created. + \param aosoa The particle AoSoA. + \param particle_created Whether to remove unused allocated space. + \param shrink_to_fit Whether to remove unused allocated space. +*/ +template +void filterEmpties( const ExecutionSpace& exec_space, + const int local_num_create, + const CreationView& particle_created, ParticleAoSoA& aosoa, + const bool shrink_to_fit ) +{ + using memory_space = typename CreationView::memory_space; + + // Determine the empty particle positions in the compaction zone. + int num_particles = aosoa.size(); + Kokkos::View empties( + Kokkos::ViewAllocateWithoutInitializing( "empties" ), + std::min( num_particles - local_num_create, local_num_create ) ); + Kokkos::parallel_scan( + "Cabana::ParticleInit::FindEmpty", + Kokkos::RangePolicy( exec_space, 0, local_num_create ), + KOKKOS_LAMBDA( const int i, int& count, const bool final_pass ) { + if ( !particle_created( i ) ) + { + if ( final_pass ) + { + empties( count ) = i; + } + ++count; + } + } ); + + // Compact the list so the it only has real particles. + Kokkos::parallel_scan( + "Cabana::ParticleInit::RemoveEmpty", + Kokkos::RangePolicy( exec_space, local_num_create, + num_particles ), + KOKKOS_LAMBDA( const int i, int& count, const bool final_pass ) { + if ( particle_created( i ) ) + { + if ( final_pass ) + { + aosoa.setTuple( empties( count ), aosoa.getTuple( i ) ); + } + ++count; + } + } ); + aosoa.resize( local_num_create ); + if ( shrink_to_fit ) + aosoa.shrinkToFit(); +} + +//---------------------------------------------------------------------------// +/*! + \brief Initialize a random number of particles in each cell given an + initialization functor. + + \tparam ParticleListType The type of particle list to initialize. + \tparam InitFunctor Initialization functor type. See the documentation below + for the create_functor parameter on the signature of this functor. + + \param exec_space Kokkos execution space. + \param create_functor A functor which populates a particle given the logical + position of a particle. This functor returns true if a particle was created + and false if it was not giving the signature: + + bool createFunctor( const double pid, const double px[3], const double pv, + typename ParticleAoSoA::tuple_type& particle ); + \param particle_list The ParticleList to populate. This will be filled with + particles and resized to a size equal to the number of particles created. + \param particles_per_cell The number of particles to sample each cell with. + \param local_grid The LocalGrid over which particles will be created. + \param shrink_to_fit Optionally remove unused allocated space after creation. + \param seed Optional random seed for generating particles. +*/ +template +void createParticles( Cabana::InitRandom, const ExecutionSpace& exec_space, + const InitFunctor& create_functor, + ParticleListType& particle_list, + const int particles_per_cell, LocalGridType& local_grid, + const bool shrink_to_fit = true, + const uint64_t seed = 123456 ) +{ + // Memory space. + using memory_space = typename ParticleListType::memory_space; + + // Create a local mesh. + auto local_mesh = Cajita::createLocalMesh( local_grid ); + + // Get the global grid. + const auto& global_grid = local_grid.globalGrid(); + + // Get the local set of owned cell indices. + auto owned_cells = + local_grid.indexSpace( Cajita::Own(), Cajita::Cell(), Cajita::Local() ); + + // Create a random number generator. + const auto local_seed = + global_grid.blockId() + ( seed % ( global_grid.blockId() + 1 ) ); + using rnd_type = Kokkos::Random_XorShift64_Pool; + rnd_type pool; + pool.init( local_seed, owned_cells.size() ); + + // Get the aosoa. + auto& aosoa = particle_list.aosoa(); + + // Allocate enough space for the case the particles consume the entire + // local grid. + int num_particles = particles_per_cell * owned_cells.size(); + aosoa.resize( num_particles ); + + // Creation status. + auto particle_created = Kokkos::View( + Kokkos::ViewAllocateWithoutInitializing( "particle_created" ), + num_particles ); + + // Initialize particles. + int local_num_create = 0; + Kokkos::parallel_reduce( + "Cajita::ParticleInit::Random", + Cajita::createExecutionPolicy( owned_cells, exec_space ), + KOKKOS_LAMBDA( const int i, const int j, const int k, + int& create_count ) { + // Compute the owned local cell id. + int i_own = i - owned_cells.min( Dim::I ); + int j_own = j - owned_cells.min( Dim::J ); + int k_own = k - owned_cells.min( Dim::K ); + int cell_id = + i_own + owned_cells.extent( Dim::I ) * + ( j_own + k_own * owned_cells.extent( Dim::J ) ); + + // Get the coordinates of the low cell node. + int low_node[3] = { i, j, k }; + double low_coords[3]; + local_mesh.coordinates( Cajita::Node(), low_node, low_coords ); + + // Get the coordinates of the high cell node. + int high_node[3] = { i + 1, j + 1, k + 1 }; + double high_coords[3]; + local_mesh.coordinates( Cajita::Node(), high_node, high_coords ); + + // Random number generator. + auto rand = pool.get_state( cell_id ); + + // Particle coordinate. + double px[3]; + + // Particle volume. + double pv = local_mesh.measure( Cajita::Cell(), low_node ) / + particles_per_cell; + + // Create particles. + for ( int p = 0; p < particles_per_cell; ++p ) + { + // Local particle id. + int pid = cell_id * particles_per_cell + p; + + // Select a random point in the cell for the particle + // location. These coordinates are logical. + for ( int d = 0; d < 3; ++d ) + { + px[d] = Kokkos::rand::draw( + rand, low_coords[d], high_coords[d] ); + } + + // Create a new particle with the given logical coordinates. + auto particle = particle_list.getParticle( pid ); + particle_created( pid ) = + create_functor( pid, px, pv, particle ); + + // If we created a new particle insert it into the list. + if ( particle_created( pid ) ) + { + particle_list.setParticle( particle, pid ); + ++create_count; + } + } + }, + local_num_create ); + + // Filter empties. + filterEmpties( exec_space, local_num_create, particle_created, aosoa, + shrink_to_fit ); +} + +/*! + \brief Initialize random particles per cell given an initialization functor. + + \tparam ParticleListType The type of particle list to initialize. + \tparam InitFunctor Initialization functor type. See the documentation below + for the create_functor parameter on the signature of this functor. + + \param tag Initialization type tag. + \param create_functor A functor which populates a particle given the logical + position of a particle. This functor returns true if a particle was created + and false if it was not giving the signature: + + bool createFunctor( const double pid, const double px[3], const double pv, + typename ParticleAoSoA::tuple_type& particle ); + \param particle_list The ParticleList to populate. This will be filled with + particles and resized to a size equal to the number of particles created. + \param particles_per_cell The number of particles to sample each cell with. + \param local_grid The LocalGrid over which particles will be created. + \param shrink_to_fit Optionally remove unused allocated space after creation. + \param seed Optional random seed for generating particles. +*/ +template +void createParticles( Cabana::InitRandom tag, const InitFunctor& create_functor, + ParticleListType& particle_list, + const int particles_per_cell, LocalGridType& local_grid, + const bool shrink_to_fit = true, + const uint64_t seed = 123456 ) +{ + using exec_space = typename ParticleListType::memory_space::execution_space; + createParticles( tag, exec_space{}, create_functor, particle_list, + particles_per_cell, local_grid, shrink_to_fit, seed ); +} + +//---------------------------------------------------------------------------// +/*! + \brief Initialize uniform particles per cell given an initialization functor. + + \tparam ParticleListType The type of particle list to initialize. + \tparam InitFunctor Initialization functor type. See the documentation below + for the create_functor parameter on the signature of this functor. + + \param exec_space Kokkos execution space. + \param create_functor A functor which populates a particle given the logical + position of a particle. This functor returns true if a particle was created + and false if it was not giving the signature: + + bool createFunctor( const double px[3], + typename ParticleAoSoA::tuple_type& particle ); + \param particle_list The ParticleList to populate. This will be filled with + particles and resized to a size equal to the number of particles created. + \param particles_per_cell_dim The number of particles to populate each cell + dimension with. + \param local_grid The LocalGrid over which particles will be created. + \param shrink_to_fit Optionally remove unused allocated space after creation. +*/ +template +void createParticles( Cabana::InitUniform, const ExecutionSpace& exec_space, + const InitFunctor& create_functor, + ParticleListType& particle_list, + const int particles_per_cell_dim, + LocalGridType& local_grid, + const bool shrink_to_fit = true ) +{ + // Memory space. + using memory_space = typename ParticleListType::memory_space; + + // Create a local mesh. + auto local_mesh = Cajita::createLocalMesh( local_grid ); + + // Get the local set of owned cell indices. + auto owned_cells = + local_grid.indexSpace( Cajita::Own(), Cajita::Cell(), Cajita::Local() ); + + // Get the aosoa. + auto& aosoa = particle_list.aosoa(); + + // Allocate enough space for particles fill the entire local grid. + int particles_per_cell = particles_per_cell_dim * particles_per_cell_dim * + particles_per_cell_dim; + int num_particles = particles_per_cell * owned_cells.size(); + aosoa.resize( num_particles ); + + // Creation status. + auto particle_created = Kokkos::View( + Kokkos::ViewAllocateWithoutInitializing( "particle_created" ), + num_particles ); + + // Initialize particles. + int local_num_create = 0; + Kokkos::parallel_reduce( + "Cajita::ParticleInit::Uniform", + Cajita::createExecutionPolicy( owned_cells, exec_space ), + KOKKOS_LAMBDA( const int i, const int j, const int k, + int& create_count ) { + // Compute the owned local cell id. + int i_own = i - owned_cells.min( Dim::I ); + int j_own = j - owned_cells.min( Dim::J ); + int k_own = k - owned_cells.min( Dim::K ); + int cell_id = + i_own + owned_cells.extent( Dim::I ) * + ( j_own + k_own * owned_cells.extent( Dim::J ) ); + + // Get the coordinates of the low cell node. + int low_node[3] = { i, j, k }; + double low_coords[3]; + local_mesh.coordinates( Cajita::Node(), low_node, low_coords ); + + // Get the coordinates of the high cell node. + int high_node[3] = { i + 1, j + 1, k + 1 }; + double high_coords[3]; + local_mesh.coordinates( Cajita::Node(), high_node, high_coords ); + + // Compute the particle spacing in each dimension. + double spacing[3] = { ( high_coords[Dim::I] - low_coords[Dim::I] ) / + particles_per_cell_dim, + ( high_coords[Dim::J] - low_coords[Dim::J] ) / + particles_per_cell_dim, + ( high_coords[Dim::K] - low_coords[Dim::K] ) / + particles_per_cell_dim }; + + // Particle coordinate. + double px[3]; + + // Particle volume. + double pv = local_mesh.measure( Cajita::Cell(), low_node ) / + particles_per_cell; + + // Create particles. + for ( int ip = 0; ip < particles_per_cell_dim; ++ip ) + for ( int jp = 0; jp < particles_per_cell_dim; ++jp ) + for ( int kp = 0; kp < particles_per_cell_dim; ++kp ) + { + // Local particle id. + int pid = cell_id * particles_per_cell + ip + + particles_per_cell_dim * + ( jp + particles_per_cell_dim * kp ); + + // Set the particle position in logical coordinates. + px[Dim::I] = 0.5 * spacing[Dim::I] + + ip * spacing[Dim::I] + low_coords[Dim::I]; + px[Dim::J] = 0.5 * spacing[Dim::J] + + jp * spacing[Dim::J] + low_coords[Dim::J]; + px[Dim::K] = 0.5 * spacing[Dim::K] + + kp * spacing[Dim::K] + low_coords[Dim::K]; + + // Create a new particle with the given logical + // coordinates. + auto particle = particle_list.getParticle( pid ); + particle_created( pid ) = + create_functor( pid, px, pv, particle ); + + // If we created a new particle insert it into the list. + if ( particle_created( pid ) ) + { + particle_list.setParticle( particle, pid ); + ++create_count; + } + } + }, + local_num_create ); + + // Filter empties. + filterEmpties( exec_space, local_num_create, particle_created, aosoa, + shrink_to_fit ); +} + +/*! + \brief Initialize random particles per cell given an initialization functor. + + \tparam ParticleListType The type of particle list to initialize. + \tparam InitFunctor Initialization functor type. See the documentation below + for the create_functor parameter on the signature of this functor. + + \param tag Initialization type tag. + \param create_functor A functor which populates a particle given the logical + position of a particle. This functor returns true if a particle was created + and false if it was not giving the signature: + + bool createFunctor( const double pid, const double px[3], const double pv, + typename ParticleAoSoA::tuple_type& particle ); + \param particle_list The ParticleList to populate. This will be filled with + particles and resized to a size equal to the number of particles created. + \param particles_per_cell The number of particles to sample each cell with. + \param local_grid The LocalGrid over which particles will be created. + \param shrink_to_fit Optionally remove unused allocated space after creation. +*/ +template +void createParticles( Cabana::InitUniform tag, + const InitFunctor& create_functor, + ParticleListType& particle_list, + const int particles_per_cell, LocalGridType& local_grid, + const bool shrink_to_fit = true ) +{ + using exec_space = typename ParticleListType::memory_space::execution_space; + createParticles( tag, exec_space{}, create_functor, particle_list, + particles_per_cell, local_grid, shrink_to_fit ); +} + +} // namespace Cajita + +#endif diff --git a/cajita/unit_test/CMakeLists.txt b/cajita/unit_test/CMakeLists.txt index 5909ae863..5dd82b18a 100644 --- a/cajita/unit_test/CMakeLists.txt +++ b/cajita/unit_test/CMakeLists.txt @@ -28,6 +28,7 @@ set(MPI_TESTS Array2d Halo3d Halo2d + ParticleInit ParticleGridDistributor2d ParticleGridDistributor3d SplineEvaluation3d diff --git a/cajita/unit_test/tstParticleInit.hpp b/cajita/unit_test/tstParticleInit.hpp new file mode 100644 index 000000000..6aa6ae83d --- /dev/null +++ b/cajita/unit_test/tstParticleInit.hpp @@ -0,0 +1,166 @@ +/**************************************************************************** + * 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 + +#include +#include +#include + +#include + +#include + +using Cajita::Dim; + +namespace Test +{ +//---------------------------------------------------------------------------// +int totalParticlesPerCell( Cabana::InitUniform, int ppc ) +{ + return ppc * ppc * ppc; +} + +int totalParticlesPerCell( Cabana::InitRandom, int ppc ) { return ppc; } + +//---------------------------------------------------------------------------// +// Field tags. +struct Foo : public Cabana::Field::Vector +{ + static std::string label() { return "foo"; } +}; + +struct Bar : public Cabana::Field::Scalar +{ + static std::string label() { return "bar"; } +}; + +//---------------------------------------------------------------------------// +template +void InitTest( InitType init_type, int ppc ) +{ + // Global bounding box. + double cell_size = 0.23; + std::array global_num_cell = { 43, 32, 39 }; + 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 ); + + std::array is_dim_periodic = { true, true, true }; + Cajita::DimBlockPartitioner<3> partitioner; + auto global_grid = Cajita::createGlobalGrid( MPI_COMM_WORLD, global_mesh, + is_dim_periodic, partitioner ); + auto local_grid = Cajita::createLocalGrid( global_grid, 0 ); + + // Make a particle list. + auto fields = Cabana::ParticleTraits(); + auto particles = + Cajita::createParticleList( "test_particles", fields ); + using plist_type = decltype( particles ); + + // Particle initialization functor. + const Kokkos::Array box = { + global_low_corner[Dim::I] + cell_size, + global_high_corner[Dim::I] - cell_size, + global_low_corner[Dim::J] + cell_size, + global_high_corner[Dim::J] - cell_size, + global_low_corner[Dim::K] + cell_size, + global_high_corner[Dim::K] - cell_size }; + auto init_func = + KOKKOS_LAMBDA( const int, const double x[3], const double v, + typename plist_type::particle_type& p ) + { + // Put particles in a box that is one cell smaller than the global + // mesh. This will give us a layer of empty cells. + if ( x[Dim::I] > box[0] && x[Dim::I] < box[1] && x[Dim::J] > box[2] && + x[Dim::J] < box[3] && x[Dim::K] > box[4] && x[Dim::K] < box[5] ) + { + for ( int d = 0; d < 3; ++d ) + get( p, Foo(), d ) = x[d]; + + get( p, Bar() ) = v; + return true; + } + else + { + return false; + } + }; + + // Initialize particles. + Cajita::createParticles( init_type, TEST_EXECSPACE(), init_func, particles, + ppc, *local_grid ); + + // Check that we made particles. + int num_p = particles.size(); + EXPECT_TRUE( num_p > 0 ); + + // Compute the global number of particles. + int global_num_particle = num_p; + MPI_Allreduce( MPI_IN_PLACE, &global_num_particle, 1, MPI_INT, MPI_SUM, + MPI_COMM_WORLD ); + int expect_num_particle = + totalParticlesPerCell( init_type, ppc ) * + ( global_grid->globalNumEntity( Cajita::Cell(), Dim::I ) - 2 ) * + ( global_grid->globalNumEntity( Cajita::Cell(), Dim::J ) - 2 ) * + ( global_grid->globalNumEntity( Cajita::Cell(), Dim::K ) - 2 ); + EXPECT_EQ( global_num_particle, expect_num_particle ); + + // Particle volume. + double cell_volume = global_mesh->cellSize( 0 ) * + global_mesh->cellSize( 1 ) * + global_mesh->cellSize( 2 ); + double volume = cell_volume / totalParticlesPerCell( init_type, ppc ); + + // Check that all particles are in the box and got initialized correctly. + auto host_particles = + Cabana::create_mirror_view_and_copy( Kokkos::HostSpace(), particles ); + for ( int p = 0; p < num_p; ++p ) + { + auto particle = host_particles.getParticle( p ); + + for ( int d = 0; d < 3; ++d ) + { + auto px = get( particle, Foo(), d ); + EXPECT_TRUE( px > box[2 * d] ); + EXPECT_TRUE( px < box[2 * d + 1] ); + } + auto pv = get( particle, Bar() ); + EXPECT_EQ( pv, volume ); + } +} + +//---------------------------------------------------------------------------// +// RUN TESTS +//---------------------------------------------------------------------------// +TEST( TEST_CATEGORY, random_init_test ) +{ + InitTest( Cabana::InitRandom(), 17 ); +} + +TEST( TEST_CATEGORY, uniform_init_test ) +{ + InitTest( Cabana::InitUniform(), 3 ); +} + +//---------------------------------------------------------------------------// + +} // end namespace Test diff --git a/core/src/Cabana_ParticleInit.hpp b/core/src/Cabana_ParticleInit.hpp index c6da3b520..c12a021e3 100644 --- a/core/src/Cabana_ParticleInit.hpp +++ b/core/src/Cabana_ParticleInit.hpp @@ -24,6 +24,64 @@ namespace Cabana { +//---------------------------------------------------------------------------// +// Initialization type tags. +//! Uniform particle initialization type tag. +struct InitUniform +{ +}; +//! Random particle initialization type tag. +struct InitRandom +{ +}; + +//---------------------------------------------------------------------------// +//! Filter out empty particles that weren't created. +template +void filterEmpties( const ExecutionSpace& exec_space, + const int local_num_create, + const CreationView& particle_created, ParticleAoSoA& aosoa ) +{ + using memory_space = typename CreationView::memory_space; + + // Determine the empty particle positions in the compaction zone. + int num_particles = aosoa.size(); + Kokkos::View empties( + Kokkos::ViewAllocateWithoutInitializing( "empties" ), + std::min( num_particles - local_num_create, local_num_create ) ); + Kokkos::parallel_scan( + "Cabana::ParticleInit::FindEmpty", + Kokkos::RangePolicy( exec_space, 0, local_num_create ), + KOKKOS_LAMBDA( const int i, int& count, const bool final_pass ) { + if ( !particle_created( i ) ) + { + if ( final_pass ) + { + empties( count ) = i; + } + ++count; + } + } ); + + // Compact the list so the it only has real particles. + Kokkos::parallel_scan( + "Cabana::ParticleInit::RemoveEmpty", + Kokkos::RangePolicy( exec_space, local_num_create, + num_particles ), + KOKKOS_LAMBDA( const int i, int& count, const bool final_pass ) { + if ( particle_created( i ) ) + { + if ( final_pass ) + { + aosoa.setTuple( empties( count ), aosoa.getTuple( i ) ); + } + ++count; + } + } ); + aosoa.resize( local_num_create ); + aosoa.shrinkToFit(); +} + /*! Generate random particles with minimum distance between neighbors. This approximates many physical scenarios, e.g. atomic simulations. Kokkos From d80e63979f788344ddfccbe57e79b71eb056ad57 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Thu, 8 Jun 2023 23:26:23 -0400 Subject: [PATCH 2/9] Rewrite core random particle init to match Cajita; Deprecate previous versions --- .../cajita/Cajita_SparseMapPerformance.cpp | 6 +- .../Cajita_SparsePartitionerPerformance.cpp | 6 +- .../core/Cabana_LinkedCellPerformance.cpp | 5 +- .../core/Cabana_NeighborArborXPerformance.cpp | 7 +- .../core/Cabana_NeighborVerletPerformance.cpp | 7 +- core/src/Cabana_ParticleInit.hpp | 604 +++++++++++++++--- core/src/Cabana_ParticleList.hpp | 17 + core/unit_test/neighbor_unit_test.hpp | 4 +- core/unit_test/tstHDF5ParticleOutput.hpp | 8 +- core/unit_test/tstParticleInit.hpp | 139 +++- 10 files changed, 661 insertions(+), 142 deletions(-) diff --git a/benchmark/cajita/Cajita_SparseMapPerformance.cpp b/benchmark/cajita/Cajita_SparseMapPerformance.cpp index 4ae216064..3e509a2bb 100644 --- a/benchmark/cajita/Cajita_SparseMapPerformance.cpp +++ b/benchmark/cajita/Cajita_SparseMapPerformance.cpp @@ -53,9 +53,9 @@ void performanceTest( std::ostream& stream, const std::string& test_prefix, positions[p] = position_type( Kokkos::ViewAllocateWithoutInitializing( "positions" ), problem_sizes[p] ); - Cabana::createRandomParticles( positions[p], problem_sizes[p], - global_low_corner[0], - global_high_corner[0] ); + Cabana::createParticles( Cabana::InitRandom(), positions[p], + problem_sizes[p], global_low_corner, + global_high_corner ); } // Number of runs in the test loops. int num_run = 10; diff --git a/benchmark/cajita/Cajita_SparsePartitionerPerformance.cpp b/benchmark/cajita/Cajita_SparsePartitionerPerformance.cpp index 41818dc9f..72ae7993f 100644 --- a/benchmark/cajita/Cajita_SparsePartitionerPerformance.cpp +++ b/benchmark/cajita/Cajita_SparsePartitionerPerformance.cpp @@ -129,9 +129,9 @@ void performanceTest( ParticleWorkloadTag, std::ostream& stream, MPI_Comm comm, positions[p] = position_type( Kokkos::ViewAllocateWithoutInitializing( "positions" ), problem_sizes[p] ); - Cabana::createRandomParticles( positions[p], problem_sizes[p], - global_low_corner[0], - global_high_corner[0] ); + Cabana::createParticles( Cabana::InitRandom(), positions[p], + problem_sizes[p], global_low_corner, + global_high_corner ); } for ( int c = 0; c < num_cells_per_dim_size; ++c ) diff --git a/benchmark/core/Cabana_LinkedCellPerformance.cpp b/benchmark/core/Cabana_LinkedCellPerformance.cpp index cd1f68bb4..64c6c9ebb 100644 --- a/benchmark/core/Cabana_LinkedCellPerformance.cpp +++ b/benchmark/core/Cabana_LinkedCellPerformance.cpp @@ -55,9 +55,12 @@ void performanceTest( std::ostream& stream, const std::string& test_prefix, // Define problem grid. x_min[p] = 0.0; x_max[p] = 1.3 * std::pow( num_p, 1.0 / 3.0 ); + double grid_min[3] = { x_min[p], x_min[p], x_min[p] }; + double grid_max[3] = { x_max[p], x_max[p], x_max[p] }; aosoas[p].resize( num_p ); auto x = Cabana::slice<0>( aosoas[p], "position" ); - Cabana::createRandomParticles( x, x.size(), x_min[p], x_max[p] ); + Cabana::createParticles( Cabana::InitRandom(), x, x.size(), grid_min, + grid_max ); } // Loop over number of ratios (neighbors per particle). diff --git a/benchmark/core/Cabana_NeighborArborXPerformance.cpp b/benchmark/core/Cabana_NeighborArborXPerformance.cpp index 9a50a9b05..36e52a55d 100644 --- a/benchmark/core/Cabana_NeighborArborXPerformance.cpp +++ b/benchmark/core/Cabana_NeighborArborXPerformance.cpp @@ -62,9 +62,12 @@ void performanceTest( std::ostream& stream, const std::string& test_prefix, // Define problem grid. x_min[p] = 0.0; x_max[p] = 1.3 * std::pow( num_p, 1.0 / 3.0 ); + double grid_min[3] = { x_min[p], x_min[p], x_min[p] }; + double grid_max[3] = { x_max[p], x_max[p], x_max[p] }; aosoas[p].resize( num_p ); auto x = Cabana::slice<0>( aosoas[p], "position" ); - Cabana::createRandomParticles( x, x.size(), x_min[p], x_max[p] ); + Cabana::createParticles( Cabana::InitRandom(), x, x.size(), grid_min, + grid_max ); if ( sort ) { @@ -74,8 +77,6 @@ void performanceTest( std::ostream& stream, const std::string& test_prefix, // in cells the size of the smallest cutoff distance. double cutoff = cutoff_ratios.front(); double sort_delta[3] = { cutoff, cutoff, cutoff }; - double grid_min[3] = { x_min[p], x_min[p], x_min[p] }; - double grid_max[3] = { x_max[p], x_max[p], x_max[p] }; auto x = Cabana::slice<0>( aosoas[p], "position" ); Cabana::LinkedCellList linked_cell_list( x, sort_delta, grid_min, grid_max ); diff --git a/benchmark/core/Cabana_NeighborVerletPerformance.cpp b/benchmark/core/Cabana_NeighborVerletPerformance.cpp index 9989b14ca..2031255c9 100644 --- a/benchmark/core/Cabana_NeighborVerletPerformance.cpp +++ b/benchmark/core/Cabana_NeighborVerletPerformance.cpp @@ -68,9 +68,12 @@ void performanceTest( std::ostream& stream, const std::string& test_prefix, // Define problem grid. x_min[p] = 0.0; x_max[p] = 1.3 * std::pow( num_p, 1.0 / 3.0 ); + double grid_min[3] = { x_min[p], x_min[p], x_min[p] }; + double grid_max[3] = { x_max[p], x_max[p], x_max[p] }; aosoas[p].resize( num_p ); auto x = Cabana::slice<0>( aosoas[p], "position" ); - Cabana::createRandomParticles( x, x.size(), x_min[p], x_max[p] ); + Cabana::createParticles( Cabana::InitRandom(), x, x.size(), grid_min, + grid_max ); if ( sort ) { @@ -80,8 +83,6 @@ void performanceTest( std::ostream& stream, const std::string& test_prefix, // in cells the size of the smallest cutoff distance. double cutoff = cutoff_ratios.front(); double sort_delta[3] = { cutoff, cutoff, cutoff }; - double grid_min[3] = { x_min[p], x_min[p], x_min[p] }; - double grid_max[3] = { x_max[p], x_max[p], x_max[p] }; Cabana::LinkedCellList linked_cell_list( x, sort_delta, grid_min, grid_max ); Cabana::permute( linked_cell_list, aosoas[p] ); diff --git a/core/src/Cabana_ParticleInit.hpp b/core/src/Cabana_ParticleInit.hpp index c12a021e3..99f1057d0 100644 --- a/core/src/Cabana_ParticleInit.hpp +++ b/core/src/Cabana_ParticleInit.hpp @@ -16,16 +16,37 @@ #ifndef CABANA_PARTICLEINIT_HPP #define CABANA_PARTICLEINIT_HPP +#include #include +#include +#include + #include #include namespace Cabana { +namespace Impl +{ +//! Copy array (std, c-array) into Kokkos::Array for potential device use. +template +auto copyArray( ArrayType corner ) +{ + using value_type = std::remove_reference_t; + Kokkos::Array kokkos_corner; + for ( std::size_t d = 0; d < 3; ++d ) + kokkos_corner[d] = corner[d]; + + return kokkos_corner; +} + +} // namespace Impl + //---------------------------------------------------------------------------// // Initialization type tags. + //! Uniform particle initialization type tag. struct InitUniform { @@ -36,69 +57,495 @@ struct InitRandom }; //---------------------------------------------------------------------------// -//! Filter out empty particles that weren't created. -template -void filterEmpties( const ExecutionSpace& exec_space, - const int local_num_create, - const CreationView& particle_created, ParticleAoSoA& aosoa ) -{ - using memory_space = typename CreationView::memory_space; - - // Determine the empty particle positions in the compaction zone. - int num_particles = aosoa.size(); - Kokkos::View empties( - Kokkos::ViewAllocateWithoutInitializing( "empties" ), - std::min( num_particles - local_num_create, local_num_create ) ); - Kokkos::parallel_scan( - "Cabana::ParticleInit::FindEmpty", - Kokkos::RangePolicy( exec_space, 0, local_num_create ), - KOKKOS_LAMBDA( const int i, int& count, const bool final_pass ) { - if ( !particle_created( i ) ) - { - if ( final_pass ) - { - empties( count ) = i; - } - ++count; - } - } ); - - // Compact the list so the it only has real particles. - Kokkos::parallel_scan( - "Cabana::ParticleInit::RemoveEmpty", - Kokkos::RangePolicy( exec_space, local_num_create, - num_particles ), - KOKKOS_LAMBDA( const int i, int& count, const bool final_pass ) { - if ( particle_created( i ) ) +// Fully random +//---------------------------------------------------------------------------// + +//---------------------------------------------------------------------------// +/*! + \brief Initialize random particles given an initialization functor. + + \param exec_space Kokkos execution space. + \param create_functor A functor which populates a particle given the logical + position of a particle. This functor returns true if a particle was created + and false if it was not giving the signature: + + bool createFunctor( const double pid, const double px[3], const double pv, + typename ParticleAoSoA::tuple_type& particle ); + \param particle_list The ParticleList to populate. This will be filled with + particles and resized to a size equal to the number of particles created. + \param num_particles The number of particles to create. + \param box_min Lower corner of volume to create particles within. Must be + device accessible if on device. + \param box_max Upper corner of volume to create particles within. Must be + device accessible if on device. + \param shrink_to_fit Optionally remove unused allocated space after creation. + \param seed Optional random seed for generating particles. + + \return Number of particles created. +*/ +template +int createParticles( + InitRandom, ExecutionSpace exec_space, const InitFunctor& create_functor, + ParticleListType& particle_list, const std::size_t num_particles, + const ArrayType box_min, const ArrayType box_max, + const bool shrink_to_fit = true, const uint64_t seed = 342343901, + typename std::enable_if::value, + int>::type* = 0 ) +{ + // Memory space. + using memory_space = typename ParticleListType::memory_space; + + using PoolType = Kokkos::Random_XorShift64_Pool; + using RandomType = Kokkos::Random_XorShift64; + PoolType pool( seed ); + + // Creation count. + auto count = Kokkos::View( "particle_count", 1 ); + + // Resize the aosoa prior to lambda capture. + auto& aosoa = particle_list.aosoa(); + aosoa.resize( num_particles ); + + // Copy corners to device accessible arrays. + auto kokkos_min = Impl::copyArray( box_min ); + auto kokkos_max = Impl::copyArray( box_max ); + + auto random_coord_op = KOKKOS_LAMBDA( const int p ) + { + // Particle coordinate. + double px[3]; + + auto gen = pool.get_state(); + auto particle = particle_list.getParticle( p ); + for ( int d = 0; d < 3; ++d ) + px[d] = Kokkos::rand::draw( gen, kokkos_min[d], + kokkos_max[d] ); + pool.free_state( gen ); + + // No volume information, so pass zero. + int create = create_functor( count( 0 ), px, 0.0, particle ); + + // If we created a new particle insert it into the list. + if ( create ) + { + auto c = Kokkos::atomic_fetch_add( &count( 0 ), 1 ); + particle_list.setParticle( particle, c ); + } + }; + + Kokkos::RangePolicy exec_policy( exec_space, 0, + num_particles ); + Kokkos::parallel_for( exec_policy, random_coord_op ); + Kokkos::fence(); + + auto count_host = + Kokkos::create_mirror_view_and_copy( Kokkos::HostSpace(), count ); + aosoa.resize( count_host( 0 ) ); + if ( shrink_to_fit ) + aosoa.shrinkToFit(); + return count_host( 0 ); +} + +//---------------------------------------------------------------------------// +/*! + \brief Initialize random particles given an initialization functor. + + \param tag Initialization type tag. + \param create_functor A functor which populates a particle given the logical + position of a particle. This functor returns true if a particle was created + and false if it was not giving the signature: + + bool createFunctor( const double pid, const double px[3], const double pv, + typename ParticleAoSoA::tuple_type& particle ); + \param particle_list The ParticleList to populate. This will be filled with + particles and resized to a size equal to the number of particles created. + \param num_particles The number of particles to create. + \param box_min Lower corner of volume to create particles within. Must be + device accessible if on device. + \param box_max Upper corner of volume to create particles within. Must be + device accessible if on device. + \param shrink_to_fit Optionally remove unused allocated space after creation. + \param seed Optional random seed for generating particles. + + \return Number of particles created. +*/ +template +int createParticles( InitRandom tag, const InitFunctor& create_functor, + ParticleListType& particle_list, + const std::size_t num_particles, const ArrayType box_min, + const ArrayType box_max, const bool shrink_to_fit = true, + const uint64_t seed = 342343901 ) +{ + using exec_space = typename ParticleListType::memory_space::execution_space; + return createParticles( tag, exec_space{}, create_functor, particle_list, + num_particles, box_min, box_max, shrink_to_fit, + seed ); +} + +/*! + \brief Initialize random particles. + + \param exec_space Kokkos execution space. + \param positions Particle positions slice. + \param num_particles The number of particles to create. + \param box_min Lower corner of volume to create particles within. Must be + device accessible if on device. + \param box_max Upper corner of volume to create particles within. Must be + device accessible if on device. + \param seed Optional random seed for generating particles. +*/ +template +void createParticles( + InitRandom, ExecutionSpace exec_space, PositionType& positions, + const std::size_t num_particles, const ArrayType box_min, + const ArrayType box_max, const uint64_t seed = 342343901, + typename std::enable_if<( is_slice::value || + Kokkos::is_view::value ), + int>::type* = 0 ) +{ + // Ensure correct space for the particles. + assert( positions.size() == num_particles ); + + // Copy corners to device accessible arrays. + auto kokkos_min = Impl::copyArray( box_min ); + auto kokkos_max = Impl::copyArray( box_max ); + + using PoolType = Kokkos::Random_XorShift64_Pool; + using RandomType = Kokkos::Random_XorShift64; + PoolType pool( seed ); + auto random_coord_op = KOKKOS_LAMBDA( const int p ) + { + auto gen = pool.get_state(); + for ( int d = 0; d < 3; ++d ) + positions( p, d ) = Kokkos::rand::draw( + gen, kokkos_min[d], kokkos_max[d] ); + pool.free_state( gen ); + }; + + Kokkos::RangePolicy exec_policy( exec_space, 0, + num_particles ); + Kokkos::parallel_for( exec_policy, random_coord_op ); + Kokkos::fence(); +} + +/*! + \brief Initialize random particles. + + \param tag Initialization type tag. + \param positions Particle positions slice. + \param num_particles The number of particles to create. + \param box_min Lower corner of volume to create particles within. Must be + device accessible if on device. + \param box_max Upper corner of volume to create particles within. Must be + device accessible if on device. + \param seed Optional random seed for generating particles. +*/ +template +void createParticles( + InitRandom tag, PositionType& positions, const std::size_t num_particles, + const ArrayType box_min, const ArrayType box_max, + const uint64_t seed = 342343901, + typename std::enable_if<( is_slice::value || + Kokkos::is_view::value ), + int>::type* = 0 ) +{ + using exec_space = typename PositionType::execution_space; + createParticles( tag, exec_space{}, positions, num_particles, box_min, + box_max, seed ); +} + +//! Generate random particles. +template +[[deprecated( "Use createParticles instead." )]] void +createRandomParticles( ExecutionSpace exec_space, PositionType& positions, + const std::size_t num_particles, const double box_min, + const double box_max ) +{ + std::array array_min = { box_min, box_min, box_min }; + std::array array_max = { box_max, box_max, box_max }; + createParticles( InitRandom{}, exec_space, positions, num_particles, + array_min, array_max ); +} + +/*! + Generate random particles. Default execution space. +*/ +template +[[deprecated( "Use createParticles instead." )]] void +createRandomParticles( PositionType& positions, const std::size_t num_particles, + const double box_min, const double box_max ) +{ + using exec_space = typename PositionType::execution_space; + createRandomParticles( exec_space{}, positions, num_particles, box_min, + box_max ); +} + +//---------------------------------------------------------------------------// +// Random with minimum separation +//---------------------------------------------------------------------------// + +/*! + \brief Initialize random particles with minimum separation. + + \param tag Initialization type tag. + \param exec_space Kokkos execution space. + \param create_functor A functor which populates a particle given the logical + position of a particle. This functor returns true if a particle was created + and false if it was not giving the signature: + + bool createFunctor( const double pid, const double px[3], const double pv, + typename ParticleAoSoA::tuple_type& particle ); + \param particle_list The ParticleList to populate. This will be filled with + particles and resized to a size equal to the number of particles created. + \param position_tag Position particle list type tag. + \param num_particles The number of particles to create. + \param min_dist Minimum separation distance between particles. Potential + particles created within this distance of an existing particle are rejected. + \param box_min Lower corner of volume to create particles within. Must be + device accessible if on device. + \param box_max Upper corner of volume to create particles within. Must be + device accessible if on device. + \param shrink_to_fit Optionally remove unused allocated space after creation. + \param seed Optional random seed for generating particles. + + \return Number of particles created. + + \note This approximates many physical scenarios, e.g. atomic simulations. +*/ +template +int createParticles( + InitRandom tag, ExecutionSpace exec_space, + const InitFunctor& create_functor, ParticleListType& particle_list, + PositionTag position_tag, const std::size_t num_particles, + const double min_dist, const ArrayType box_min, const ArrayType box_max, + const bool shrink_to_fit = true, const uint64_t seed = 342343901, + typename std::enable_if::value, + int>::type* = 0 ) +{ + double min_dist_sqr = min_dist * min_dist; + + // Resize the aosoa prior to lambda capture. + auto& aosoa = particle_list.aosoa(); + aosoa.resize( num_particles ); + + auto positions = particle_list.slice( position_tag ); + + // Create the functor which ignores particles within the radius of another. + auto min_distance_op = + KOKKOS_LAMBDA( const int id, const double px[3], const double, + typename ParticleListType::particle_type& particle ) + { + // Ensure this particle is not within the minimum distance of any other + // existing particle. + for ( int n = 0; n < id; n++ ) + { + double dx = positions( n, 0 ) - px[0]; + double dy = positions( n, 1 ) - px[1]; + double dz = positions( n, 2 ) - px[2]; + double dist = dx * dx + dy * dy + dz * dz; + if ( dist < min_dist_sqr ) + return false; + } + + bool create = create_functor( id, px, 0.0, particle ); + return create; + }; + + // Pass the functor to the general case. + return createParticles( tag, exec_space, min_distance_op, particle_list, + num_particles, box_min, box_max, shrink_to_fit, + seed ); +} + +/*! + \brief Initialize random particles with minimum separation. + + \param tag Initialization type tag. + \param create_functor A functor which populates a particle given the logical + position of a particle. This functor returns true if a particle was created + and false if it was not giving the signature: + + bool createFunctor( const double pid, const double px[3], const double pv, + typename ParticleAoSoA::tuple_type& particle ); + \param particle_list The ParticleList to populate. This will be filled with + particles and resized to a size equal to the number of particles created. + \param position_tag Position particle list type tag. + \param num_particles The number of particles to create. + \param min_dist Minimum separation distance between particles. Potential + particles created within this distance of an existing particle are rejected. + \param box_min Lower corner of volume to create particles within. Must be + device accessible if on device. + \param box_max Upper corner of volume to create particles within. Must be + device accessible if on device. + \param shrink_to_fit Optionally remove unused allocated space after creation. + \param seed Optional random seed for generating particles. + + \return Number of particles created. + + \note This approximates many physical scenarios, e.g. atomic simulations. +*/ +template +int createParticles( InitRandom tag, const InitFunctor& create_functor, + ParticleListType& particle_list, PositionTag position_tag, + const std::size_t num_particles, const double min_dist, + const ArrayType box_min, const ArrayType box_max, + const bool shrink_to_fit = true, + const uint64_t seed = 342343901 ) +{ + using exec_space = typename ParticleListType::memory_space::execution_space; + return createParticles( tag, exec_space{}, create_functor, particle_list, + position_tag, num_particles, min_dist, box_min, + box_max, shrink_to_fit, seed ); +} + +/*! + \brief Initialize random particles with minimum separation. + + \param exec_space Kokkos execution space. + \param positions Particle positions slice. + \param num_particles The number of particles to create. + \param min_dist Minimum separation distance between particles. Potential + particles created within this distance of an existing particle are rejected. + \param box_min Lower corner of volume to create particles within. Must be + device accessible if on device. + \param box_max Upper corner of volume to create particles within. Must be + device accessible if on device. + \param seed Optional random seed for generating particles. + + \return Number of particles created. + + \note This approximates many physical scenarios, e.g. atomic simulations. +*/ +template +int createParticles( + InitRandom, ExecutionSpace exec_space, PositionType& positions, + const std::size_t num_particles, const double min_dist, + const ArrayType box_min, const ArrayType box_max, + const uint64_t seed = 342343901, + typename std::enable_if::value, int>::type* = 0 ) +{ + double min_dist_sqr = min_dist * min_dist; + + using PoolType = Kokkos::Random_XorShift64_Pool; + using RandomType = Kokkos::Random_XorShift64; + PoolType pool( seed ); + + // Creation count. + using memory_space = typename PositionType::memory_space; + auto count = Kokkos::View( "particle_count", 1 ); + + auto random_coord_op = KOKKOS_LAMBDA( const int p ) + { + auto gen = pool.get_state(); + + double px[3]; + for ( int d = 0; d < 3; ++d ) + px[d] = Kokkos::rand::draw( gen, box_min[d], + box_max[d] ); + pool.free_state( gen ); + + for ( int n = 0; n < p; n++ ) + { + double dx = positions( n, 0 ) - px[0]; + double dy = positions( n, 1 ) - px[1]; + double dz = positions( n, 2 ) - px[2]; + double dist = dx * dx + dy * dy + dz * dz; + if ( dist < min_dist_sqr ) { - if ( final_pass ) - { - aosoa.setTuple( empties( count ), aosoa.getTuple( i ) ); - } - ++count; + return; } - } ); - aosoa.resize( local_num_create ); - aosoa.shrinkToFit(); + } + // Only increment if no neighbors within min distance. + auto c = Kokkos::atomic_fetch_add( &count( 0 ), 1 ); + for ( int d = 0; d < 3; ++d ) + positions( c, d ) = px[d]; + }; + + Kokkos::RangePolicy exec_policy( exec_space, 0, + num_particles ); + Kokkos::parallel_for( exec_policy, random_coord_op ); + Kokkos::fence(); + + auto count_host = + Kokkos::create_mirror_view_and_copy( Kokkos::HostSpace(), count ); + return count_host( 0 ); } /*! - Generate random particles with minimum distance between neighbors. This - approximates many physical scenarios, e.g. atomic simulations. Kokkos - device version. + \brief Initialize random particles with minimum separation. + + \param tag Initialization type tag. + \param positions Particle positions slice. + \param num_particles The number of particles to create. + \param min_dist Minimum separation distance between particles. Potential + particles created within this distance of an existing particle are rejected. + \param box_min Lower corner of volume to create particles within. Must be + device accessible if on device. + \param box_max Upper corner of volume to create particles within. Must be + device accessible if on device. + \param seed Optional random seed for generating particles. + + \return Number of particles created. + + \note This approximates many physical scenarios, e.g. atomic simulations. +*/ +template +int createParticles( + InitRandom tag, PositionType& positions, const std::size_t num_particles, + const double min_dist, const ArrayType box_max, const ArrayType box_min, + const uint64_t seed = 342343901, + typename std::enable_if::value, int>::type* = 0 ) +{ + using exec_space = typename PositionType::execution_space; + return createParticles( tag, exec_space{}, positions, num_particles, + min_dist, box_min, box_max, seed ); +} + +/*! + \brief Generate random particles with minimum distance between neighbors. + \note This approximates many physical scenarios, e.g. atomic simulations. + + \note This version will continue sampling until the number of selected + particles is created. This can be extremely slow based on the requested box + size and minimum separation distance. +*/ +template +[[deprecated( "Use createParticles instead." )]] void +createRandomParticlesMinDistance( PositionType& positions, + const std::size_t num_particles, + const double box_min, const double box_max, + const double min_dist ) +{ + using exec_space = typename PositionType::execution_space; + createRandomParticlesMinDistance( exec_space{}, positions, num_particles, + box_min, box_max, min_dist ); +} + +/*! + \brief Generate random particles with minimum distance between neighbors. + \note This approximates many physical scenarios, e.g. atomic simulations. + + This version will continue sampling until the number of selected particles is + created. This can be extremely slow based on the requested box size and + minimum separation distance. */ template -void createRandomParticlesMinDistance( ExecutionSpace, PositionType& positions, - const std::size_t num_particles, - const double box_min, - const double box_max, - const double min_dist ) +[[deprecated( "Use createParticles instead." )]] void +createRandomParticlesMinDistance( ExecutionSpace exec_space, + PositionType& positions, + const std::size_t num_particles, + const double box_min, const double box_max, + const double min_dist ) { double min_dist_sqr = min_dist * min_dist; using PoolType = Kokkos::Random_XorShift64_Pool; using RandomType = Kokkos::Random_XorShift64; PoolType pool( 342343901 ); + auto random_coord_op = KOKKOS_LAMBDA( const int p ) { auto gen = pool.get_state(); @@ -129,67 +576,16 @@ void createRandomParticlesMinDistance( ExecutionSpace, PositionType& positions, break; } } - pool.free_state( gen ); } - }; - Kokkos::RangePolicy exec_policy( 0, num_particles ); - Kokkos::parallel_for( exec_policy, random_coord_op ); - Kokkos::fence(); -} - -/*! - Generate random particles with minimum distance between neighbors. This - approximates many physical scenarios, e.g. atomic simulations. Kokkos - device version with default execution space. -*/ -template -void createRandomParticlesMinDistance( PositionType& positions, - const std::size_t num_particles, - const double box_min, - const double box_max, - const double min_dist ) -{ - using exec_space = typename PositionType::execution_space; - createRandomParticlesMinDistance( exec_space{}, positions, num_particles, - box_min, box_max, min_dist ); -} - -//! Generate random particles. Kokkos device version. -template -void createRandomParticles( ExecutionSpace, PositionType& positions, - const std::size_t num_particles, - const double box_min, const double box_max ) -{ - using PoolType = Kokkos::Random_XorShift64_Pool; - using RandomType = Kokkos::Random_XorShift64; - PoolType pool( 342343901 ); - auto random_coord_op = KOKKOS_LAMBDA( const int p ) - { - auto gen = pool.get_state(); - for ( int d = 0; d < 3; ++d ) - positions( p, d ) = - Kokkos::rand::draw( gen, box_min, box_max ); pool.free_state( gen ); }; - Kokkos::RangePolicy exec_policy( 0, num_particles ); + + Kokkos::RangePolicy exec_policy( exec_space, 0, + num_particles ); Kokkos::parallel_for( exec_policy, random_coord_op ); Kokkos::fence(); } -/*! - Generate random particles. Kokkos device version with default execution - space. -*/ -template -void createRandomParticles( PositionType& positions, - const std::size_t num_particles, - const double box_min, const double box_max ) -{ - using exec_space = typename PositionType::execution_space; - createRandomParticles( exec_space{}, positions, num_particles, box_min, - box_max ); -} - } // namespace Cabana #endif diff --git a/core/src/Cabana_ParticleList.hpp b/core/src/Cabana_ParticleList.hpp index 7f3269d3f..7d22f4911 100644 --- a/core/src/Cabana_ParticleList.hpp +++ b/core/src/Cabana_ParticleList.hpp @@ -262,6 +262,23 @@ class ParticleList aosoa_type _aosoa; }; +template +struct is_particle_list_impl : public std::false_type +{ +}; + +template +struct is_particle_list_impl> + : public std::true_type +{ +}; + +//! ParticleList static type checker. +template +struct is_particle_list + : public is_particle_list_impl::type>::type +{ +}; //---------------------------------------------------------------------------// /*! \brief ParticleList creation function. diff --git a/core/unit_test/neighbor_unit_test.hpp b/core/unit_test/neighbor_unit_test.hpp index 28ff84257..30cf55875 100644 --- a/core/unit_test/neighbor_unit_test.hpp +++ b/core/unit_test/neighbor_unit_test.hpp @@ -965,8 +965,8 @@ struct NeighborListTestData auto positions = Cabana::slice<0>( aosoa ); #endif - Cabana::createRandomParticles( positions, positions.size(), box_min, - box_max ); + Cabana::createParticles( Cabana::InitRandom(), positions, + positions.size(), grid_min, grid_max ); #ifdef KOKKOS_ENABLE_OPENMPTARGET // FIXME_OPENMPTARGET Cabana::deep_copy( aosoa, aosoa_copy ); diff --git a/core/unit_test/tstHDF5ParticleOutput.hpp b/core/unit_test/tstHDF5ParticleOutput.hpp index 0b97dc2d6..5e3e8c4aa 100644 --- a/core/unit_test/tstHDF5ParticleOutput.hpp +++ b/core/unit_test/tstHDF5ParticleOutput.hpp @@ -72,8 +72,8 @@ void checkMatrix( SliceType1 write, SliceType2 read ) //---------------------------------------------------------------------------// void writeReadTest() { - double low_corner = -2.8; - double high_corner = 1.2; + std::array low_corner = { -2.8, 1.4, -10.4 }; + std::array high_corner = { 1.2, 7.5, -7.9 }; // Allocate particle properties. int num_particle = 100; @@ -88,8 +88,8 @@ void writeReadTest() auto ids = Cabana::slice<3>( aosoa, "ids" ); // Create random particles. - Cabana::createRandomParticles( coords, num_particle, low_corner, - high_corner ); + Cabana::createParticles( Cabana::InitRandom(), coords, num_particle, + low_corner, high_corner ); // Set other particle properties. auto aosoa_mirror = diff --git a/core/unit_test/tstParticleInit.hpp b/core/unit_test/tstParticleInit.hpp index 4e2a13008..f4dd778c6 100644 --- a/core/unit_test/tstParticleInit.hpp +++ b/core/unit_test/tstParticleInit.hpp @@ -17,9 +17,23 @@ #include +namespace Test +{ + +struct Foo : Cabana::Field::Scalar +{ + static std::string label() { return "foo"; } +}; + +struct Bar : Cabana::Field::Matrix +{ + static std::string label() { return "bar"; } +}; + template -void checkRandomParticles( const int num_particle, const double box_min, - const double box_max, +void checkRandomParticles( const int num_particle, + const Kokkos::Array box_min, + const Kokkos::Array box_max, const PositionType host_positions ) { // Check that we got as many particles as we should have. @@ -29,8 +43,8 @@ void checkRandomParticles( const int num_particle, const double box_min, for ( int p = 0; p < num_particle; ++p ) for ( int d = 0; d < 3; ++d ) { - EXPECT_GE( host_positions( p, d ), box_min ); - EXPECT_LE( host_positions( p, d ), box_max ); + EXPECT_GE( host_positions( p, d ), box_min[d] ); + EXPECT_LE( host_positions( p, d ), box_max[d] ); } } @@ -50,11 +64,11 @@ void checkRandomDistances( const int min_distance, double diff = host_positions( i, d ) - host_positions( j, d ); dsqr += diff * diff; } - EXPECT_GE( dsqr, min_distance ); + EXPECT_GE( dsqr, min_distance * min_distance ); } } -void testRandomCreationMinDistance() +void testRandomCreationSliceMinDistance() { int num_particle = 200; Cabana::AoSoA, TEST_MEMSPACE> aosoa( @@ -62,29 +76,36 @@ void testRandomCreationMinDistance() auto positions = Cabana::slice<0>( aosoa ); double min_dist = 0.47; - double box_min = -9.5; - double box_max = 7.6; - Cabana::createRandomParticlesMinDistance( positions, positions.size(), - box_min, box_max, min_dist ); + Kokkos::Array box_min = { -9.5, -4.7, 0.5 }; + Kokkos::Array box_max = { 7.6, -1.5, 5.5 }; + int created = + Cabana::createParticles( Cabana::InitRandom(), positions, + positions.size(), min_dist, box_min, box_max ); + aosoa.resize( created ); auto host_aosoa = Cabana::create_mirror_view_and_copy( Kokkos::HostSpace(), aosoa ); auto host_positions = Cabana::slice<0>( host_aosoa ); - checkRandomParticles( num_particle, box_min, box_max, host_positions ); + // All particles may not have been created in this case (some skipped by the + // minimum distance criterion). + EXPECT_LE( host_positions.size(), num_particle ); + EXPECT_GE( host_positions.size(), 0 ); + checkRandomParticles( host_positions.size(), box_min, box_max, + host_positions ); checkRandomDistances( min_dist, host_positions ); } -void testRandomCreation() +void testRandomCreationSlice() { int num_particle = 200; Cabana::AoSoA, TEST_MEMSPACE> aosoa( "random", num_particle ); auto positions = Cabana::slice<0>( aosoa ); - double box_min = -9.5; - double box_max = 7.6; - Cabana::createRandomParticles( positions, positions.size(), box_min, - box_max ); + Kokkos::Array box_min = { -9.5, -4.7, 0.5 }; + Kokkos::Array box_max = { 7.6, -1.5, 5.5 }; + Cabana::createParticles( Cabana::InitRandom(), positions, positions.size(), + box_min, box_max ); auto host_aosoa = Cabana::create_mirror_view_and_copy( Kokkos::HostSpace(), aosoa ); auto host_positions = Cabana::slice<0>( host_aosoa ); @@ -92,8 +113,88 @@ void testRandomCreation() checkRandomParticles( num_particle, box_min, box_max, host_positions ); } -TEST( TEST_CATEGORY, random_particle_creation_test ) +void testRandomCreationParticleListMinDistance() { - testRandomCreationMinDistance(); - testRandomCreation(); + int num_particle = 200; + using plist_type = + Cabana::ParticleList, Foo, + Bar>; + plist_type particle_list( "random_particles" ); + + double min_dist = 0.47; + Kokkos::Array box_min = { -9.5, -4.7, 0.5 }; + Kokkos::Array box_max = { 7.6, -1.5, 5.5 }; + // Create all particles. + auto init_func = + KOKKOS_LAMBDA( const int, const double x[3], const double, + typename plist_type::particle_type& particle ) + { + // Set positions. + for ( int d = 0; d < 3; ++d ) + Cabana::get( particle, Cabana::Field::Position<3>(), d ) = x[d]; + + return true; + }; + Cabana::createParticles( Cabana::InitRandom(), init_func, particle_list, + Cabana::Field::Position<3>(), num_particle, + min_dist, box_min, box_max ); + + auto host_particle_list = Cabana::create_mirror_view_and_copy( + Kokkos::HostSpace(), particle_list ); + auto host_positions = + host_particle_list.slice( Cabana::Field::Position<3>() ); + + // All particles may not have been created in this case (some skipped by the + // minimum distance criterion). + EXPECT_LE( host_positions.size(), num_particle ); + EXPECT_GE( host_positions.size(), 0 ); + checkRandomParticles( host_particle_list.size(), box_min, box_max, + host_positions ); + checkRandomDistances( min_dist, host_positions ); } + +void testRandomCreationParticleList() +{ + int num_particle = 200; + using plist_type = + Cabana::ParticleList, Foo, + Bar>; + plist_type particle_list( "random_particles" ); + + Kokkos::Array box_min = { -9.5, -4.7, 0.5 }; + Kokkos::Array box_max = { 7.6, -1.5, 5.5 }; + // Create all particles. + auto init_func = + KOKKOS_LAMBDA( const int, const double x[3], const double, + typename plist_type::particle_type& particle ) + { + for ( int d = 0; d < 3; ++d ) + Cabana::get( particle, Cabana::Field::Position<3>(), d ) = x[d]; + + return true; + }; + auto created = + Cabana::createParticles( Cabana::InitRandom(), init_func, particle_list, + num_particle, box_min, box_max ); + EXPECT_EQ( num_particle, created ); + + auto host_particle_list = Cabana::create_mirror_view_and_copy( + Kokkos::HostSpace(), particle_list ); + auto host_positions = + host_particle_list.slice( Cabana::Field::Position<3>() ); + + checkRandomParticles( num_particle, box_min, box_max, host_positions ); +} + +TEST( TEST_CATEGORY, random_particle_creation_slice_test ) +{ + testRandomCreationSliceMinDistance(); + testRandomCreationSlice(); +} +TEST( TEST_CATEGORY, random_particle_creation_particlelist_test ) +{ + testRandomCreationParticleListMinDistance(); + testRandomCreationParticleList(); +} + +} // namespace Test \ No newline at end of file From 6611e05a25123623fa349e693f5e96bea3dbb6d1 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 26 Jun 2023 14:50:14 -0400 Subject: [PATCH 3/9] Add Cajita particle init slice variant --- cajita/src/Cajita_ParticleInit.hpp | 313 +++++++++++++++++++++++---- cajita/src/Cajita_ParticleList.hpp | 18 ++ cajita/unit_test/tstParticleInit.hpp | 72 +++++- 3 files changed, 353 insertions(+), 50 deletions(-) diff --git a/cajita/src/Cajita_ParticleInit.hpp b/cajita/src/Cajita_ParticleInit.hpp index b817535da..f92093e21 100644 --- a/cajita/src/Cajita_ParticleInit.hpp +++ b/cajita/src/Cajita_ParticleInit.hpp @@ -35,6 +35,8 @@ #include #include +#include + namespace Cajita { @@ -100,10 +102,6 @@ void filterEmpties( const ExecutionSpace& exec_space, \brief Initialize a random number of particles in each cell given an initialization functor. - \tparam ParticleListType The type of particle list to initialize. - \tparam InitFunctor Initialization functor type. See the documentation below - for the create_functor parameter on the signature of this functor. - \param exec_space Kokkos execution space. \param create_functor A functor which populates a particle given the logical position of a particle. This functor returns true if a particle was created @@ -120,12 +118,13 @@ void filterEmpties( const ExecutionSpace& exec_space, */ template -void createParticles( Cabana::InitRandom, const ExecutionSpace& exec_space, - const InitFunctor& create_functor, - ParticleListType& particle_list, - const int particles_per_cell, LocalGridType& local_grid, - const bool shrink_to_fit = true, - const uint64_t seed = 123456 ) +void createParticles( + Cabana::InitRandom, const ExecutionSpace& exec_space, + const InitFunctor& create_functor, ParticleListType& particle_list, + const int particles_per_cell, LocalGridType& local_grid, + const bool shrink_to_fit = true, const uint64_t seed = 123456, + typename std::enable_if::value, + int>::type* = 0 ) { // Memory space. using memory_space = typename ParticleListType::memory_space; @@ -162,9 +161,8 @@ void createParticles( Cabana::InitRandom, const ExecutionSpace& exec_space, // Initialize particles. int local_num_create = 0; - Kokkos::parallel_reduce( - "Cajita::ParticleInit::Random", - Cajita::createExecutionPolicy( owned_cells, exec_space ), + Cajita::grid_parallel_reduce( + "Cajita::ParticleInit::Random", exec_space, owned_cells, KOKKOS_LAMBDA( const int i, const int j, const int k, int& create_count ) { // Compute the owned local cell id. @@ -251,11 +249,13 @@ void createParticles( Cabana::InitRandom, const ExecutionSpace& exec_space, \param seed Optional random seed for generating particles. */ template -void createParticles( Cabana::InitRandom tag, const InitFunctor& create_functor, - ParticleListType& particle_list, - const int particles_per_cell, LocalGridType& local_grid, - const bool shrink_to_fit = true, - const uint64_t seed = 123456 ) +void createParticles( + Cabana::InitRandom tag, const InitFunctor& create_functor, + ParticleListType& particle_list, const int particles_per_cell, + LocalGridType& local_grid, const bool shrink_to_fit = true, + const uint64_t seed = 123456, + typename std::enable_if::value, + int>::type* = 0 ) { using exec_space = typename ParticleListType::memory_space::execution_space; createParticles( tag, exec_space{}, create_functor, particle_list, @@ -264,21 +264,131 @@ void createParticles( Cabana::InitRandom tag, const InitFunctor& create_functor, //---------------------------------------------------------------------------// /*! - \brief Initialize uniform particles per cell given an initialization functor. + \brief Initialize a random number of particles in each cell. + + \param exec_space Kokkos execution space. + \param positions Particle positions slice. This should be already the size of + the number of grid cells times particles_per_cell.s + \param particles_per_cell The number of particles to sample each cell with. + \param local_grid The LocalGrid over which particles will be created. + \param seed Optional random seed for generating particles. +*/ +template +void createParticles( + Cabana::InitRandom, const ExecutionSpace& exec_space, + PositionType& positions, const int particles_per_cell, + LocalGridType& local_grid, const uint64_t seed = 123456, + typename std::enable_if<( Cabana::is_slice::value || + Kokkos::is_view::value ), + int>::type* = 0 ) +{ + // Create a local mesh. + auto local_mesh = Cajita::createLocalMesh( local_grid ); + + // Get the global grid. + const auto& global_grid = local_grid.globalGrid(); + + // Get the local set of owned cell indices. + auto owned_cells = + local_grid.indexSpace( Cajita::Own(), Cajita::Cell(), Cajita::Local() ); + + // Create a random number generator. + const auto local_seed = + global_grid.blockId() + ( seed % ( global_grid.blockId() + 1 ) ); + using rnd_type = Kokkos::Random_XorShift64_Pool; + rnd_type pool; + pool.init( local_seed, owned_cells.size() ); + + // Ensure correct space for the particles. + assert( positions.size() == static_cast( + particles_per_cell * owned_cells.size() ) ); + + // Initialize particles. + Cajita::grid_parallel_for( + "Cajita::ParticleInit::Random", exec_space, owned_cells, + KOKKOS_LAMBDA( const int i, const int j, const int k ) { + // Compute the owned local cell id. + int i_own = i - owned_cells.min( Dim::I ); + int j_own = j - owned_cells.min( Dim::J ); + int k_own = k - owned_cells.min( Dim::K ); + int cell_id = + i_own + owned_cells.extent( Dim::I ) * + ( j_own + k_own * owned_cells.extent( Dim::J ) ); + + // Get the coordinates of the low cell node. + int low_node[3] = { i, j, k }; + double low_coords[3]; + local_mesh.coordinates( Cajita::Node(), low_node, low_coords ); + + // Get the coordinates of the high cell node. + int high_node[3] = { i + 1, j + 1, k + 1 }; + double high_coords[3]; + local_mesh.coordinates( Cajita::Node(), high_node, high_coords ); + + // Random number generator. + auto rand = pool.get_state( cell_id ); + + // Create particles. + for ( int p = 0; p < particles_per_cell; ++p ) + { + // Local particle id. + int pid = cell_id * particles_per_cell + p; + + // Select a random point in the cell for the particle + // location. These coordinates are logical. + for ( int d = 0; d < 3; ++d ) + { + positions( pid, d ) = + Kokkos::rand::draw( + rand, low_coords[d], high_coords[d] ); + } + } + } ); +} + +/*! + \brief Initialize a random number of particles in each cell. + + \param tag Initialization type tag. + \param positions Particle positions slice. This should be already the size of + the number of grid cells times particles_per_cell.s + \param particles_per_cell The number of particles to sample each cell with. + \param local_grid The LocalGrid over which particles will be created. + \param seed Optional random seed for generating particles. +*/ +template +void createParticles( + Cabana::InitRandom tag, PositionType& positions, + const int particles_per_cell, LocalGridType& local_grid, + const uint64_t seed = 123456, + typename std::enable_if<( Cabana::is_slice::value || + Kokkos::is_view::value ), + int>::type* = 0 ) +{ + using exec_space = typename PositionType::execution_space; + createParticles( tag, exec_space{}, positions, particles_per_cell, + local_grid, seed ); +} + +//---------------------------------------------------------------------------// +/*! + \brief Initialize uniform particles per cell given an initialization + functor. \tparam ParticleListType The type of particle list to initialize. - \tparam InitFunctor Initialization functor type. See the documentation below - for the create_functor parameter on the signature of this functor. + \tparam InitFunctor Initialization functor type. See the documentation + below for the create_functor parameter on the signature of this functor. \param exec_space Kokkos execution space. - \param create_functor A functor which populates a particle given the logical - position of a particle. This functor returns true if a particle was created - and false if it was not giving the signature: + \param create_functor A functor which populates a particle given the + logical position of a particle. This functor returns true if a particle + was created and false if it was not giving the signature: bool createFunctor( const double px[3], typename ParticleAoSoA::tuple_type& particle ); - \param particle_list The ParticleList to populate. This will be filled with - particles and resized to a size equal to the number of particles created. + \param particle_list The ParticleList to populate. This will be filled + with particles and resized to a size equal to the number of particles + created. \param particles_per_cell_dim The number of particles to populate each cell dimension with. \param local_grid The LocalGrid over which particles will be created. @@ -286,12 +396,13 @@ void createParticles( Cabana::InitRandom tag, const InitFunctor& create_functor, */ template -void createParticles( Cabana::InitUniform, const ExecutionSpace& exec_space, - const InitFunctor& create_functor, - ParticleListType& particle_list, - const int particles_per_cell_dim, - LocalGridType& local_grid, - const bool shrink_to_fit = true ) +void createParticles( + Cabana::InitUniform, const ExecutionSpace& exec_space, + const InitFunctor& create_functor, ParticleListType& particle_list, + const int particles_per_cell_dim, LocalGridType& local_grid, + const bool shrink_to_fit = true, + typename std::enable_if::value, + int>::type* = 0 ) { // Memory space. using memory_space = typename ParticleListType::memory_space; @@ -319,9 +430,8 @@ void createParticles( Cabana::InitUniform, const ExecutionSpace& exec_space, // Initialize particles. int local_num_create = 0; - Kokkos::parallel_reduce( - "Cajita::ParticleInit::Uniform", - Cajita::createExecutionPolicy( owned_cells, exec_space ), + Cajita::grid_parallel_reduce( + "Cajita::ParticleInit::Uniform", exec_space, owned_cells, KOKKOS_LAMBDA( const int i, const int j, const int k, int& create_count ) { // Compute the owned local cell id. @@ -397,11 +507,7 @@ void createParticles( Cabana::InitUniform, const ExecutionSpace& exec_space, } /*! - \brief Initialize random particles per cell given an initialization functor. - - \tparam ParticleListType The type of particle list to initialize. - \tparam InitFunctor Initialization functor type. See the documentation below - for the create_functor parameter on the signature of this functor. + \brief Initialize uniform particles per cell given an initialization functor. \param tag Initialization type tag. \param create_functor A functor which populates a particle given the logical @@ -412,22 +518,135 @@ void createParticles( Cabana::InitUniform, const ExecutionSpace& exec_space, typename ParticleAoSoA::tuple_type& particle ); \param particle_list The ParticleList to populate. This will be filled with particles and resized to a size equal to the number of particles created. - \param particles_per_cell The number of particles to sample each cell with. + \param particles_per_cell_dim The number of particles to populate each cell + dimension with. \param local_grid The LocalGrid over which particles will be created. \param shrink_to_fit Optionally remove unused allocated space after creation. */ template -void createParticles( Cabana::InitUniform tag, - const InitFunctor& create_functor, - ParticleListType& particle_list, - const int particles_per_cell, LocalGridType& local_grid, - const bool shrink_to_fit = true ) +void createParticles( + Cabana::InitUniform tag, const InitFunctor& create_functor, + ParticleListType& particle_list, const int particles_per_cell_dim, + LocalGridType& local_grid, const bool shrink_to_fit = true, + typename std::enable_if::value, + int>::type* = 0 ) { using exec_space = typename ParticleListType::memory_space::execution_space; createParticles( tag, exec_space{}, create_functor, particle_list, - particles_per_cell, local_grid, shrink_to_fit ); + particles_per_cell_dim, local_grid, shrink_to_fit ); +} + +//---------------------------------------------------------------------------// +/*! + \brief Initialize a uniform number of particles in each cell. + + \param exec_space Kokkos execution space. + \param positions Particle positions slice. This should be already the size of + the number of grid cells times particles_per_cell.s + \param particles_per_cell_dim The number of particles to populate each cell + dimension with. + \param local_grid The LocalGrid over which particles will be created. +*/ +template +void createParticles( + Cabana::InitUniform, const ExecutionSpace& exec_space, + PositionType& positions, const int particles_per_cell_dim, + LocalGridType& local_grid, + typename std::enable_if<( Cabana::is_slice::value || + Kokkos::is_view::value ), + int>::type* = 0 ) +{ + // Create a local mesh. + auto local_mesh = Cajita::createLocalMesh( local_grid ); + + // Get the local set of owned cell indices. + auto owned_cells = + local_grid.indexSpace( Cajita::Own(), Cajita::Cell(), Cajita::Local() ); + + int particles_per_cell = particles_per_cell_dim * particles_per_cell_dim * + particles_per_cell_dim; + + // Ensure correct space for the particles. + assert( positions.size() == static_cast( + particles_per_cell * owned_cells.size() ) ); + + // Initialize particles. + Cajita::grid_parallel_for( + "Cajita::ParticleInit::Uniform", exec_space, owned_cells, + KOKKOS_LAMBDA( const int i, const int j, const int k ) { + // Compute the owned local cell id. + int i_own = i - owned_cells.min( Dim::I ); + int j_own = j - owned_cells.min( Dim::J ); + int k_own = k - owned_cells.min( Dim::K ); + int cell_id = + i_own + owned_cells.extent( Dim::I ) * + ( j_own + k_own * owned_cells.extent( Dim::J ) ); + + // Get the coordinates of the low cell node. + int low_node[3] = { i, j, k }; + double low_coords[3]; + local_mesh.coordinates( Cajita::Node(), low_node, low_coords ); + + // Get the coordinates of the high cell node. + int high_node[3] = { i + 1, j + 1, k + 1 }; + double high_coords[3]; + local_mesh.coordinates( Cajita::Node(), high_node, high_coords ); + + // Compute the particle spacing in each dimension. + double spacing[3] = { ( high_coords[Dim::I] - low_coords[Dim::I] ) / + particles_per_cell_dim, + ( high_coords[Dim::J] - low_coords[Dim::J] ) / + particles_per_cell_dim, + ( high_coords[Dim::K] - low_coords[Dim::K] ) / + particles_per_cell_dim }; + + // Create particles. + for ( int ip = 0; ip < particles_per_cell_dim; ++ip ) + for ( int jp = 0; jp < particles_per_cell_dim; ++jp ) + for ( int kp = 0; kp < particles_per_cell_dim; ++kp ) + { + // Local particle id. + int pid = cell_id * particles_per_cell + ip + + particles_per_cell_dim * + ( jp + particles_per_cell_dim * kp ); + + // Set the particle position in logical coordinates. + positions( pid, 0 ) = 0.5 * spacing[Dim::I] + + ip * spacing[Dim::I] + + low_coords[Dim::I]; + positions( pid, 1 ) = 0.5 * spacing[Dim::J] + + jp * spacing[Dim::J] + + low_coords[Dim::J]; + positions( pid, 2 ) = 0.5 * spacing[Dim::K] + + kp * spacing[Dim::K] + + low_coords[Dim::K]; + } + } ); } +//---------------------------------------------------------------------------// +/*! + \brief Initialize a uniform number of particles in each cell. + + \param tag Initialization type tag. + \param positions Particle positions slice. This should be already the size of + the number of grid cells times particles_per_cell.s + \param particles_per_cell_dim The number of particles to populate each cell + dimension with. + \param local_grid The LocalGrid over which particles will be created. +*/ +template +void createParticles( + Cabana::InitUniform tag, PositionType& positions, + const int particles_per_cell_dim, LocalGridType& local_grid, + typename std::enable_if<( Cabana::is_slice::value || + Kokkos::is_view::value ), + int>::type* = 0 ) +{ + using exec_space = typename PositionType::execution_space; + createParticles( tag, exec_space{}, positions, particles_per_cell_dim, + local_grid ); +} } // namespace Cajita #endif diff --git a/cajita/src/Cajita_ParticleList.hpp b/cajita/src/Cajita_ParticleList.hpp index 16955a457..a8dafb243 100644 --- a/cajita/src/Cajita_ParticleList.hpp +++ b/cajita/src/Cajita_ParticleList.hpp @@ -97,6 +97,24 @@ class ParticleList : public Cabana::ParticleList using base::_aosoa; }; +template +struct is_particle_list_impl : public std::false_type +{ +}; + +template +struct is_particle_list_impl> + : public std::true_type +{ +}; + +//! ParticleList static type checker. +template +struct is_particle_list + : public is_particle_list_impl::type>::type +{ +}; + //---------------------------------------------------------------------------// /*! \brief ParticleList creation function. diff --git a/cajita/unit_test/tstParticleInit.hpp b/cajita/unit_test/tstParticleInit.hpp index 6aa6ae83d..982cf77a5 100644 --- a/cajita/unit_test/tstParticleInit.hpp +++ b/cajita/unit_test/tstParticleInit.hpp @@ -51,7 +51,7 @@ struct Bar : public Cabana::Field::Scalar //---------------------------------------------------------------------------// template -void InitTest( InitType init_type, int ppc ) +void initParticleListTest( InitType init_type, int ppc ) { // Global bounding box. double cell_size = 0.23; @@ -148,17 +148,83 @@ void InitTest( InitType init_type, int ppc ) } } +//---------------------------------------------------------------------------// +template +void initSliceTest( InitType init_type, int ppc ) +{ + // Global bounding box. + double cell_size = 0.23; + std::array global_num_cell = { 43, 32, 39 }; + 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 ); + + std::array is_dim_periodic = { true, true, true }; + Cajita::DimBlockPartitioner<3> partitioner; + auto global_grid = Cajita::createGlobalGrid( MPI_COMM_WORLD, global_mesh, + is_dim_periodic, partitioner ); + auto local_grid = Cajita::createLocalGrid( global_grid, 0 ); + auto owned_cells = local_grid->indexSpace( Cajita::Own(), Cajita::Cell(), + Cajita::Local() ); + + int num_particle = + owned_cells.size() * totalParticlesPerCell( init_type, ppc ); + Cabana::AoSoA, TEST_MEMSPACE> aosoa( + "random", num_particle ); + auto positions = Cabana::slice<0>( aosoa ); + + // Particle initialization functor. + const Kokkos::Array box = { + global_low_corner[Dim::I], global_high_corner[Dim::I], + global_low_corner[Dim::J], global_high_corner[Dim::J], + global_low_corner[Dim::K], global_high_corner[Dim::K] }; + + // Initialize all particles. + Cajita::createParticles( init_type, TEST_EXECSPACE(), positions, ppc, + *local_grid ); + + // Check that we created all particles. + int global_num_particle = positions.size(); + MPI_Allreduce( MPI_IN_PLACE, &global_num_particle, 1, MPI_INT, MPI_SUM, + MPI_COMM_WORLD ); + int expect_num_particle = + totalParticlesPerCell( init_type, ppc ) * + global_grid->globalNumEntity( Cajita::Cell(), Dim::I ) * + global_grid->globalNumEntity( Cajita::Cell(), Dim::J ) * + global_grid->globalNumEntity( Cajita::Cell(), Dim::K ); + EXPECT_EQ( global_num_particle, expect_num_particle ); + + // Check that all particles are in the box. + auto host_aosoa = + Cabana::create_mirror_view_and_copy( Kokkos::HostSpace(), aosoa ); + auto host_positions = Cabana::slice<0>( host_aosoa ); + for ( std::size_t p = 0; p < host_positions.size(); ++p ) + { + for ( std::size_t d = 0; d < 3; ++d ) + { + EXPECT_TRUE( host_positions( p, d ) > box[2 * d] ); + EXPECT_TRUE( host_positions( p, d ) < box[2 * d + 1] ); + } + } +} + //---------------------------------------------------------------------------// // RUN TESTS //---------------------------------------------------------------------------// TEST( TEST_CATEGORY, random_init_test ) { - InitTest( Cabana::InitRandom(), 17 ); + initParticleListTest( Cabana::InitRandom(), 17 ); + initSliceTest( Cabana::InitRandom(), 17 ); } TEST( TEST_CATEGORY, uniform_init_test ) { - InitTest( Cabana::InitUniform(), 3 ); + initParticleListTest( Cabana::InitUniform(), 3 ); + initSliceTest( Cabana::InitUniform(), 3 ); } //---------------------------------------------------------------------------// From 30ef034e627eedf69ccdad8d6596edcd08ac5cbd Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Tue, 27 Jun 2023 12:51:41 -0400 Subject: [PATCH 4/9] fixup: Cajita particle init return num created --- cajita/src/Cajita_ParticleInit.hpp | 27 +++++++++++---------------- 1 file changed, 11 insertions(+), 16 deletions(-) diff --git a/cajita/src/Cajita_ParticleInit.hpp b/cajita/src/Cajita_ParticleInit.hpp index f92093e21..42ce74844 100644 --- a/cajita/src/Cajita_ParticleInit.hpp +++ b/cajita/src/Cajita_ParticleInit.hpp @@ -118,7 +118,7 @@ void filterEmpties( const ExecutionSpace& exec_space, */ template -void createParticles( +int createParticles( Cabana::InitRandom, const ExecutionSpace& exec_space, const InitFunctor& create_functor, ParticleListType& particle_list, const int particles_per_cell, LocalGridType& local_grid, @@ -225,15 +225,12 @@ void createParticles( // Filter empties. filterEmpties( exec_space, local_num_create, particle_created, aosoa, shrink_to_fit ); + return local_num_create; } /*! \brief Initialize random particles per cell given an initialization functor. - \tparam ParticleListType The type of particle list to initialize. - \tparam InitFunctor Initialization functor type. See the documentation below - for the create_functor parameter on the signature of this functor. - \param tag Initialization type tag. \param create_functor A functor which populates a particle given the logical position of a particle. This functor returns true if a particle was created @@ -249,7 +246,7 @@ void createParticles( \param seed Optional random seed for generating particles. */ template -void createParticles( +int createParticles( Cabana::InitRandom tag, const InitFunctor& create_functor, ParticleListType& particle_list, const int particles_per_cell, LocalGridType& local_grid, const bool shrink_to_fit = true, @@ -258,8 +255,9 @@ void createParticles( int>::type* = 0 ) { using exec_space = typename ParticleListType::memory_space::execution_space; - createParticles( tag, exec_space{}, create_functor, particle_list, - particles_per_cell, local_grid, shrink_to_fit, seed ); + return createParticles( tag, exec_space{}, create_functor, particle_list, + particles_per_cell, local_grid, shrink_to_fit, + seed ); } //---------------------------------------------------------------------------// @@ -375,10 +373,6 @@ void createParticles( \brief Initialize uniform particles per cell given an initialization functor. - \tparam ParticleListType The type of particle list to initialize. - \tparam InitFunctor Initialization functor type. See the documentation - below for the create_functor parameter on the signature of this functor. - \param exec_space Kokkos execution space. \param create_functor A functor which populates a particle given the logical position of a particle. This functor returns true if a particle @@ -396,7 +390,7 @@ void createParticles( */ template -void createParticles( +int createParticles( Cabana::InitUniform, const ExecutionSpace& exec_space, const InitFunctor& create_functor, ParticleListType& particle_list, const int particles_per_cell_dim, LocalGridType& local_grid, @@ -504,6 +498,7 @@ void createParticles( // Filter empties. filterEmpties( exec_space, local_num_create, particle_created, aosoa, shrink_to_fit ); + return local_num_create; } /*! @@ -524,7 +519,7 @@ void createParticles( \param shrink_to_fit Optionally remove unused allocated space after creation. */ template -void createParticles( +int createParticles( Cabana::InitUniform tag, const InitFunctor& create_functor, ParticleListType& particle_list, const int particles_per_cell_dim, LocalGridType& local_grid, const bool shrink_to_fit = true, @@ -532,8 +527,8 @@ void createParticles( int>::type* = 0 ) { using exec_space = typename ParticleListType::memory_space::execution_space; - createParticles( tag, exec_space{}, create_functor, particle_list, - particles_per_cell_dim, local_grid, shrink_to_fit ); + return createParticles( tag, exec_space{}, create_functor, particle_list, + particles_per_cell_dim, local_grid, shrink_to_fit ); } //---------------------------------------------------------------------------// From 91b36ddc0f72174b7e89fba37eebfd3d396e3559 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 30 Jun 2023 16:14:11 -0400 Subject: [PATCH 5/9] Remove new slice min distance variant; size will be wrong and would need filtering --- core/src/Cabana_ParticleInit.hpp | 103 ----------------------------- core/unit_test/tstParticleInit.hpp | 28 -------- 2 files changed, 131 deletions(-) diff --git a/core/src/Cabana_ParticleInit.hpp b/core/src/Cabana_ParticleInit.hpp index 99f1057d0..a89797227 100644 --- a/core/src/Cabana_ParticleInit.hpp +++ b/core/src/Cabana_ParticleInit.hpp @@ -401,109 +401,6 @@ int createParticles( InitRandom tag, const InitFunctor& create_functor, box_max, shrink_to_fit, seed ); } -/*! - \brief Initialize random particles with minimum separation. - - \param exec_space Kokkos execution space. - \param positions Particle positions slice. - \param num_particles The number of particles to create. - \param min_dist Minimum separation distance between particles. Potential - particles created within this distance of an existing particle are rejected. - \param box_min Lower corner of volume to create particles within. Must be - device accessible if on device. - \param box_max Upper corner of volume to create particles within. Must be - device accessible if on device. - \param seed Optional random seed for generating particles. - - \return Number of particles created. - - \note This approximates many physical scenarios, e.g. atomic simulations. -*/ -template -int createParticles( - InitRandom, ExecutionSpace exec_space, PositionType& positions, - const std::size_t num_particles, const double min_dist, - const ArrayType box_min, const ArrayType box_max, - const uint64_t seed = 342343901, - typename std::enable_if::value, int>::type* = 0 ) -{ - double min_dist_sqr = min_dist * min_dist; - - using PoolType = Kokkos::Random_XorShift64_Pool; - using RandomType = Kokkos::Random_XorShift64; - PoolType pool( seed ); - - // Creation count. - using memory_space = typename PositionType::memory_space; - auto count = Kokkos::View( "particle_count", 1 ); - - auto random_coord_op = KOKKOS_LAMBDA( const int p ) - { - auto gen = pool.get_state(); - - double px[3]; - for ( int d = 0; d < 3; ++d ) - px[d] = Kokkos::rand::draw( gen, box_min[d], - box_max[d] ); - pool.free_state( gen ); - - for ( int n = 0; n < p; n++ ) - { - double dx = positions( n, 0 ) - px[0]; - double dy = positions( n, 1 ) - px[1]; - double dz = positions( n, 2 ) - px[2]; - double dist = dx * dx + dy * dy + dz * dz; - if ( dist < min_dist_sqr ) - { - return; - } - } - // Only increment if no neighbors within min distance. - auto c = Kokkos::atomic_fetch_add( &count( 0 ), 1 ); - for ( int d = 0; d < 3; ++d ) - positions( c, d ) = px[d]; - }; - - Kokkos::RangePolicy exec_policy( exec_space, 0, - num_particles ); - Kokkos::parallel_for( exec_policy, random_coord_op ); - Kokkos::fence(); - - auto count_host = - Kokkos::create_mirror_view_and_copy( Kokkos::HostSpace(), count ); - return count_host( 0 ); -} - -/*! - \brief Initialize random particles with minimum separation. - - \param tag Initialization type tag. - \param positions Particle positions slice. - \param num_particles The number of particles to create. - \param min_dist Minimum separation distance between particles. Potential - particles created within this distance of an existing particle are rejected. - \param box_min Lower corner of volume to create particles within. Must be - device accessible if on device. - \param box_max Upper corner of volume to create particles within. Must be - device accessible if on device. - \param seed Optional random seed for generating particles. - - \return Number of particles created. - - \note This approximates many physical scenarios, e.g. atomic simulations. -*/ -template -int createParticles( - InitRandom tag, PositionType& positions, const std::size_t num_particles, - const double min_dist, const ArrayType box_max, const ArrayType box_min, - const uint64_t seed = 342343901, - typename std::enable_if::value, int>::type* = 0 ) -{ - using exec_space = typename PositionType::execution_space; - return createParticles( tag, exec_space{}, positions, num_particles, - min_dist, box_min, box_max, seed ); -} - /*! \brief Generate random particles with minimum distance between neighbors. \note This approximates many physical scenarios, e.g. atomic simulations. diff --git a/core/unit_test/tstParticleInit.hpp b/core/unit_test/tstParticleInit.hpp index f4dd778c6..a4fd552b9 100644 --- a/core/unit_test/tstParticleInit.hpp +++ b/core/unit_test/tstParticleInit.hpp @@ -68,33 +68,6 @@ void checkRandomDistances( const int min_distance, } } -void testRandomCreationSliceMinDistance() -{ - int num_particle = 200; - Cabana::AoSoA, TEST_MEMSPACE> aosoa( - "random", num_particle ); - auto positions = Cabana::slice<0>( aosoa ); - - double min_dist = 0.47; - Kokkos::Array box_min = { -9.5, -4.7, 0.5 }; - Kokkos::Array box_max = { 7.6, -1.5, 5.5 }; - int created = - Cabana::createParticles( Cabana::InitRandom(), positions, - positions.size(), min_dist, box_min, box_max ); - aosoa.resize( created ); - auto host_aosoa = - Cabana::create_mirror_view_and_copy( Kokkos::HostSpace(), aosoa ); - auto host_positions = Cabana::slice<0>( host_aosoa ); - - // All particles may not have been created in this case (some skipped by the - // minimum distance criterion). - EXPECT_LE( host_positions.size(), num_particle ); - EXPECT_GE( host_positions.size(), 0 ); - checkRandomParticles( host_positions.size(), box_min, box_max, - host_positions ); - checkRandomDistances( min_dist, host_positions ); -} - void testRandomCreationSlice() { int num_particle = 200; @@ -188,7 +161,6 @@ void testRandomCreationParticleList() TEST( TEST_CATEGORY, random_particle_creation_slice_test ) { - testRandomCreationSliceMinDistance(); testRandomCreationSlice(); } TEST( TEST_CATEGORY, random_particle_creation_particlelist_test ) From d79896060639f22f103d33e3a45177b2856404f0 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 30 Jun 2023 16:44:46 -0400 Subject: [PATCH 6/9] Separate checks for slice and View particle sizes --- core/src/Cabana_ParticleInit.hpp | 4 ++-- core/src/Cabana_Slice.hpp | 20 ++++++++++++++++++++ 2 files changed, 22 insertions(+), 2 deletions(-) diff --git a/core/src/Cabana_ParticleInit.hpp b/core/src/Cabana_ParticleInit.hpp index a89797227..f597586c9 100644 --- a/core/src/Cabana_ParticleInit.hpp +++ b/core/src/Cabana_ParticleInit.hpp @@ -204,8 +204,8 @@ void createParticles( Kokkos::is_view::value ), int>::type* = 0 ) { - // Ensure correct space for the particles. - assert( positions.size() == num_particles ); + // Ensure correct space for the particles (View or Slice). + checkSize( positions, num_particles ); // Copy corners to device accessible arrays. auto kokkos_min = Impl::copyArray( box_min ); diff --git a/core/src/Cabana_Slice.hpp b/core/src/Cabana_Slice.hpp index 899f5d668..587f80620 100644 --- a/core/src/Cabana_Slice.hpp +++ b/core/src/Cabana_Slice.hpp @@ -981,6 +981,26 @@ void copyViewToSlice( ViewType& view, const SliceType& slice, copyViewToSlice( exec_space{}, view, slice, begin, end ); } +//! Check slice size (differs from Kokkos View). +template +void checkSize( + [[maybe_unused]] SliceType slice, [[maybe_unused]] const std::size_t size, + typename std::enable_if::value, int>::type* = 0 ) +{ + // Allow unused for release builds. + assert( slice.size() == size ); +} + +//! Check View size (differs from Slice). +template +void checkSize( + [[maybe_unused]] ViewType view, [[maybe_unused]] const std::size_t size, + typename std::enable_if::value, int>::type* = 0 ) +{ + // Allow unused for release builds. + assert( view.extent( 0 ) == size ); +} + //---------------------------------------------------------------------------// } // end namespace Cabana From 005fe2b153f23fdbdea16fdf07f5b33afd6faaa5 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 2 Aug 2023 09:40:04 -0400 Subject: [PATCH 7/9] fixup: create aosoa with size in particle list deep copy --- core/src/Cabana_DeepCopy.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/core/src/Cabana_DeepCopy.hpp b/core/src/Cabana_DeepCopy.hpp index 96a69364a..993947f16 100644 --- a/core/src/Cabana_DeepCopy.hpp +++ b/core/src/Cabana_DeepCopy.hpp @@ -438,7 +438,8 @@ auto create_mirror_view_and_copy( // Create an AoSoA in the new memory space. using src_plist_type = ParticleList; using member_types = typename src_plist_type::member_types; - AoSoA aosoa_dst( aosoa_src.label() ); + AoSoA aosoa_dst( aosoa_src.label(), + aosoa_src.size() ); // Copy data to new AoAoA. deep_copy( aosoa_dst, aosoa_src ); From d101d22e8e4b045733c8c2a5de06609d0557219e Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 2 Aug 2023 11:27:53 -0400 Subject: [PATCH 8/9] fixup: particle list AoSoA constructor --- core/src/Cabana_ParticleList.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/core/src/Cabana_ParticleList.hpp b/core/src/Cabana_ParticleList.hpp index 7d22f4911..e1099500c 100644 --- a/core/src/Cabana_ParticleList.hpp +++ b/core/src/Cabana_ParticleList.hpp @@ -207,8 +207,7 @@ class ParticleList } //! Constructor from existing AoSoA. - template - ParticleList( const AoSoAType aosoa ) + ParticleList( const aosoa_type& aosoa ) : _aosoa( aosoa ) { } From a3763e5f0674952cd427479a38eaaa1650440df9 Mon Sep 17 00:00:00 2001 From: Samuel Reeve Date: Mon, 14 Aug 2023 19:26:02 +0000 Subject: [PATCH 9/9] Convert from fill with filter empty to atomic grid create --- cajita/src/Cajita_ParticleInit.hpp | 129 +++++++-------------------- cajita/unit_test/tstParticleInit.hpp | 9 +- 2 files changed, 38 insertions(+), 100 deletions(-) diff --git a/cajita/src/Cajita_ParticleInit.hpp b/cajita/src/Cajita_ParticleInit.hpp index 42ce74844..4867b3a20 100644 --- a/cajita/src/Cajita_ParticleInit.hpp +++ b/cajita/src/Cajita_ParticleInit.hpp @@ -40,63 +40,6 @@ namespace Cajita { -//---------------------------------------------------------------------------// -/*! - \brief Filter out empty particles that weren't created. - - \param exec_space Kokkos execution space. - \param local_num_create The number of particles created. - \param aosoa The particle AoSoA. - \param particle_created Whether to remove unused allocated space. - \param shrink_to_fit Whether to remove unused allocated space. -*/ -template -void filterEmpties( const ExecutionSpace& exec_space, - const int local_num_create, - const CreationView& particle_created, ParticleAoSoA& aosoa, - const bool shrink_to_fit ) -{ - using memory_space = typename CreationView::memory_space; - - // Determine the empty particle positions in the compaction zone. - int num_particles = aosoa.size(); - Kokkos::View empties( - Kokkos::ViewAllocateWithoutInitializing( "empties" ), - std::min( num_particles - local_num_create, local_num_create ) ); - Kokkos::parallel_scan( - "Cabana::ParticleInit::FindEmpty", - Kokkos::RangePolicy( exec_space, 0, local_num_create ), - KOKKOS_LAMBDA( const int i, int& count, const bool final_pass ) { - if ( !particle_created( i ) ) - { - if ( final_pass ) - { - empties( count ) = i; - } - ++count; - } - } ); - - // Compact the list so the it only has real particles. - Kokkos::parallel_scan( - "Cabana::ParticleInit::RemoveEmpty", - Kokkos::RangePolicy( exec_space, local_num_create, - num_particles ), - KOKKOS_LAMBDA( const int i, int& count, const bool final_pass ) { - if ( particle_created( i ) ) - { - if ( final_pass ) - { - aosoa.setTuple( empties( count ), aosoa.getTuple( i ) ); - } - ++count; - } - } ); - aosoa.resize( local_num_create ); - if ( shrink_to_fit ) - aosoa.shrinkToFit(); -} - //---------------------------------------------------------------------------// /*! \brief Initialize a random number of particles in each cell given an @@ -154,17 +97,13 @@ int createParticles( int num_particles = particles_per_cell * owned_cells.size(); aosoa.resize( num_particles ); - // Creation status. - auto particle_created = Kokkos::View( - Kokkos::ViewAllocateWithoutInitializing( "particle_created" ), - num_particles ); + // Creation count. + auto count = Kokkos::View( "particle_count", 1 ); // Initialize particles. - int local_num_create = 0; - Cajita::grid_parallel_reduce( + Cajita::grid_parallel_for( "Cajita::ParticleInit::Random", exec_space, owned_cells, - KOKKOS_LAMBDA( const int i, const int j, const int k, - int& create_count ) { + KOKKOS_LAMBDA( const int i, const int j, const int k ) { // Compute the owned local cell id. int i_own = i - owned_cells.min( Dim::I ); int j_own = j - owned_cells.min( Dim::J ); @@ -209,23 +148,24 @@ int createParticles( // Create a new particle with the given logical coordinates. auto particle = particle_list.getParticle( pid ); - particle_created( pid ) = - create_functor( pid, px, pv, particle ); + bool created = create_functor( pid, px, pv, particle ); // If we created a new particle insert it into the list. - if ( particle_created( pid ) ) + if ( created ) { - particle_list.setParticle( particle, pid ); - ++create_count; + auto c = Kokkos::atomic_fetch_add( &count( 0 ), 1 ); + particle_list.setParticle( particle, c ); } } - }, - local_num_create ); + } ); + Kokkos::fence(); - // Filter empties. - filterEmpties( exec_space, local_num_create, particle_created, aosoa, - shrink_to_fit ); - return local_num_create; + auto count_host = + Kokkos::create_mirror_view_and_copy( Kokkos::HostSpace(), count ); + aosoa.resize( count_host( 0 ) ); + if ( shrink_to_fit ) + aosoa.shrinkToFit(); + return count_host( 0 ); } /*! @@ -417,17 +357,13 @@ int createParticles( int num_particles = particles_per_cell * owned_cells.size(); aosoa.resize( num_particles ); - // Creation status. - auto particle_created = Kokkos::View( - Kokkos::ViewAllocateWithoutInitializing( "particle_created" ), - num_particles ); + // Creation count. + auto count = Kokkos::View( "particle_count", 1 ); // Initialize particles. - int local_num_create = 0; - Cajita::grid_parallel_reduce( + Cajita::grid_parallel_for( "Cajita::ParticleInit::Uniform", exec_space, owned_cells, - KOKKOS_LAMBDA( const int i, const int j, const int k, - int& create_count ) { + KOKKOS_LAMBDA( const int i, const int j, const int k ) { // Compute the owned local cell id. int i_own = i - owned_cells.min( Dim::I ); int j_own = j - owned_cells.min( Dim::J ); @@ -482,23 +418,26 @@ int createParticles( // Create a new particle with the given logical // coordinates. auto particle = particle_list.getParticle( pid ); - particle_created( pid ) = - create_functor( pid, px, pv, particle ); + bool created = create_functor( pid, px, pv, particle ); // If we created a new particle insert it into the list. - if ( particle_created( pid ) ) + if ( created ) { - particle_list.setParticle( particle, pid ); - ++create_count; + // TODO: replace atomic with fill (use pid directly) + // and filter of empty list entries. + auto c = Kokkos::atomic_fetch_add( &count( 0 ), 1 ); + particle_list.setParticle( particle, c ); } } - }, - local_num_create ); + } ); + Kokkos::fence(); - // Filter empties. - filterEmpties( exec_space, local_num_create, particle_created, aosoa, - shrink_to_fit ); - return local_num_create; + auto count_host = + Kokkos::create_mirror_view_and_copy( Kokkos::HostSpace(), count ); + aosoa.resize( count_host( 0 ) ); + if ( shrink_to_fit ) + aosoa.shrinkToFit(); + return count_host( 0 ); } /*! diff --git a/cajita/unit_test/tstParticleInit.hpp b/cajita/unit_test/tstParticleInit.hpp index 982cf77a5..dc45aa695 100644 --- a/cajita/unit_test/tstParticleInit.hpp +++ b/cajita/unit_test/tstParticleInit.hpp @@ -106,15 +106,14 @@ void initParticleListTest( InitType init_type, int ppc ) }; // Initialize particles. - Cajita::createParticles( init_type, TEST_EXECSPACE(), init_func, particles, - ppc, *local_grid ); + int num_p = Cajita::createParticles( init_type, TEST_EXECSPACE(), init_func, + particles, ppc, *local_grid ); // Check that we made particles. - int num_p = particles.size(); EXPECT_TRUE( num_p > 0 ); // Compute the global number of particles. - int global_num_particle = num_p; + int global_num_particle = particles.size(); MPI_Allreduce( MPI_IN_PLACE, &global_num_particle, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD ); int expect_num_particle = @@ -144,7 +143,7 @@ void initParticleListTest( InitType init_type, int ppc ) EXPECT_TRUE( px < box[2 * d + 1] ); } auto pv = get( particle, Bar() ); - EXPECT_EQ( pv, volume ); + EXPECT_DOUBLE_EQ( pv, volume ); } }