Cabana 0.8.0-dev
 
Loading...
Searching...
No Matches
Cabana_ParticleInit.hpp
Go to the documentation of this file.
1/****************************************************************************
2 * Copyright (c) 2018-2023 by the Cabana authors *
3 * All rights reserved. *
4 * *
5 * This file is part of the Cabana library. Cabana is distributed under a *
6 * BSD 3-clause license. For the licensing terms see the LICENSE file in *
7 * the top-level directory. *
8 * *
9 * SPDX-License-Identifier: BSD-3-Clause *
10 ****************************************************************************/
11
16#ifndef CABANA_PARTICLEINIT_HPP
17#define CABANA_PARTICLEINIT_HPP
18
19#include <Kokkos_Core.hpp>
20#include <Kokkos_Random.hpp>
21
22#include <Cabana_ParticleList.hpp>
23#include <Cabana_Slice.hpp>
24#include <Cabana_Utils.hpp>
25
26#include <random>
27#include <type_traits>
28
29namespace Cabana
30{
31
32//---------------------------------------------------------------------------//
33// Initialization type tags.
34
36struct InitUniform
37{
38};
39
40struct InitRandom
41{
42};
43
44//---------------------------------------------------------------------------//
45// Fully random
46//---------------------------------------------------------------------------//
47
48//---------------------------------------------------------------------------//
71template <class ExecutionSpace, class InitFunctor, class ParticleListType,
72 template <class, std::size_t, class...> class ArrayType, class Scalar,
73 std::size_t NumSpaceDim, class... Args>
74int createParticles(
75 InitRandom, ExecutionSpace exec_space, const InitFunctor& create_functor,
76 ParticleListType& particle_list, const std::size_t num_particles,
77#ifdef KOKKOS_ENABLE_DEPRECATED_CODE_4
78 const ArrayType<Scalar, NumSpaceDim, Args...> box_min,
79 const ArrayType<Scalar, NumSpaceDim, Args...> box_max,
80#else
81 const ArrayType<Scalar, NumSpaceDim> box_min,
82 const ArrayType<Scalar, NumSpaceDim> box_max,
83#endif
84 const std::size_t previous_num_particles = 0,
85 const bool shrink_to_fit = true, const uint64_t seed = 342343901,
86 typename std::enable_if<is_particle_list<ParticleListType>::value,
87 int>::type* = 0 )
88{
89 // Memory space.
90 using memory_space = typename ParticleListType::memory_space;
91
92 using PoolType = Kokkos::Random_XorShift64_Pool<ExecutionSpace>;
93 using RandomType = Kokkos::Random_XorShift64<ExecutionSpace>;
94 PoolType pool( seed );
95
96 // Resize the aosoa prior to lambda capture.
97 auto& aosoa = particle_list.aosoa();
98 aosoa.resize( previous_num_particles + num_particles );
99
100 // Creation count.
101 auto count = Kokkos::View<int*, memory_space>( "particle_count", 1 );
102 Kokkos::deep_copy( count, previous_num_particles );
103
104 // Copy corners to device accessible arrays.
105 auto kokkos_min = copyArray( box_min );
106 auto kokkos_max = copyArray( box_max );
107
108 auto random_coord_op = KOKKOS_LAMBDA( const int p )
109 {
110 // Particle coordinate.
111 double px[NumSpaceDim];
112
113 auto gen = pool.get_state();
114 auto particle = particle_list.getParticle( p );
115 for ( std::size_t d = 0; d < NumSpaceDim; ++d )
116 px[d] = Kokkos::rand<RandomType, double>::draw( gen, kokkos_min[d],
117 kokkos_max[d] );
118 pool.free_state( gen );
119
120 // No volume information, so pass zero.
121 int create = create_functor( count( 0 ), px, 0.0, particle );
122
123 // If we created a new particle insert it into the list.
124 if ( create )
125 {
126 auto c = Kokkos::atomic_fetch_add( &count( 0 ), 1 );
127 particle_list.setParticle( particle, c );
128 }
129 };
130
131 Kokkos::RangePolicy<ExecutionSpace> exec_policy(
132 exec_space, previous_num_particles,
133 previous_num_particles + num_particles );
134 Kokkos::parallel_for( "Cabana::createParticles", exec_policy,
135 random_coord_op );
136 Kokkos::fence();
137
138 auto count_host =
139 Kokkos::create_mirror_view_and_copy( Kokkos::HostSpace(), count );
140 aosoa.resize( count_host( 0 ) );
141 if ( shrink_to_fit )
142 aosoa.shrinkToFit();
143 return count_host( 0 );
144}
145
146//---------------------------------------------------------------------------//
169template <class InitFunctor, class ParticleListType,
170 template <class, std::size_t, class...> class ArrayType, class Scalar,
171 std::size_t NumSpaceDim, class... Args>
172int createParticles( InitRandom tag, const InitFunctor& create_functor,
173 ParticleListType& particle_list,
174 const std::size_t num_particles,
175#ifdef KOKKOS_ENABLE_DEPRECATED_CODE_4
176 const ArrayType<Scalar, NumSpaceDim, Args...> box_min,
177 const ArrayType<Scalar, NumSpaceDim, Args...> box_max,
178#else
179 const ArrayType<Scalar, NumSpaceDim> box_min,
180 const ArrayType<Scalar, NumSpaceDim> box_max,
181#endif
182 const std::size_t previous_num_particles = 0,
183 const bool shrink_to_fit = true,
184 const uint64_t seed = 342343901 )
185{
186 using exec_space = typename ParticleListType::memory_space::execution_space;
187 return createParticles( tag, exec_space{}, create_functor, particle_list,
188 num_particles, box_min, box_max,
189 previous_num_particles, shrink_to_fit, seed );
190}
191
204template <class ExecutionSpace, class PositionType,
205 template <class, std::size_t, class...> class ArrayType, class Scalar,
206 std::size_t NumSpaceDim, class... Args>
207void createParticles(
208 InitRandom, ExecutionSpace exec_space, PositionType& positions,
209 const std::size_t num_particles,
210#ifdef KOKKOS_ENABLE_DEPRECATED_CODE_4
211 const ArrayType<Scalar, NumSpaceDim, Args...> box_min,
212 const ArrayType<Scalar, NumSpaceDim, Args...> box_max,
213#else
214 const ArrayType<Scalar, NumSpaceDim> box_min,
215 const ArrayType<Scalar, NumSpaceDim> box_max,
216#endif
217 const std::size_t previous_num_particles = 0,
218 const uint64_t seed = 342343901,
219 typename std::enable_if<( is_slice<PositionType>::value ||
220 Kokkos::is_view<PositionType>::value ),
221 int>::type* = 0 )
222{
223 // Ensure correct space for the particles (View or Slice).
224 checkSize( positions, num_particles + previous_num_particles );
225
226 // Copy corners to device accessible arrays.
227 auto kokkos_min = copyArray( box_min );
228 auto kokkos_max = copyArray( box_max );
229
230 using PoolType = Kokkos::Random_XorShift64_Pool<ExecutionSpace>;
231 using RandomType = Kokkos::Random_XorShift64<ExecutionSpace>;
232 PoolType pool( seed );
233 auto random_coord_op = KOKKOS_LAMBDA( const int p )
234 {
235 auto gen = pool.get_state();
236 for ( std::size_t d = 0; d < NumSpaceDim; ++d )
237 positions( p, d ) = Kokkos::rand<RandomType, double>::draw(
238 gen, kokkos_min[d], kokkos_max[d] );
239 pool.free_state( gen );
240 };
241
242 Kokkos::RangePolicy<ExecutionSpace> exec_policy(
243 exec_space, previous_num_particles,
244 previous_num_particles + num_particles );
245 Kokkos::parallel_for( "Cabana::createParticles", exec_policy,
246 random_coord_op );
247 Kokkos::fence();
248}
249
262template <class PositionType,
263 template <class, std::size_t, class...> class ArrayType, class Scalar,
264 std::size_t NumSpaceDim, class... Args>
265void createParticles(
266 InitRandom tag, PositionType& positions, const std::size_t num_particles,
267#ifdef KOKKOS_ENABLE_DEPRECATED_CODE_4
268 const ArrayType<Scalar, NumSpaceDim, Args...> box_min,
269 const ArrayType<Scalar, NumSpaceDim, Args...> box_max,
270#else
271 const ArrayType<Scalar, NumSpaceDim> box_min,
272 const ArrayType<Scalar, NumSpaceDim> box_max,
273#endif
274 const std::size_t previous_num_particles = 0,
275 const uint64_t seed = 342343901,
276 typename std::enable_if<( is_slice<PositionType>::value ||
277 Kokkos::is_view<PositionType>::value ),
278 int>::type* = 0 )
279{
280 using exec_space = typename PositionType::execution_space;
281 createParticles( tag, exec_space{}, positions, num_particles, box_min,
282 box_max, previous_num_particles, seed );
283}
284
298template <class ExecutionSpace, class PositionType, class Scalar>
299void createParticles( InitRandom tag, ExecutionSpace exec_space,
300 PositionType& positions, const std::size_t num_particles,
301 const Scalar box_min[3], const Scalar box_max[3],
302 const std::size_t previous_num_particles = 0,
303 const uint64_t seed = 342343901 )
304{
305 auto kokkos_min = copyArray<Scalar, 3>( box_min );
306 auto kokkos_max = copyArray<Scalar, 3>( box_max );
307 createParticles( tag, exec_space, positions, num_particles, kokkos_min,
308 kokkos_max, previous_num_particles, seed );
309}
310
323template <class PositionType, class Scalar>
324void createParticles( InitRandom tag, PositionType& positions,
325 const std::size_t num_particles, const Scalar box_min[3],
326 const Scalar box_max[3],
327 const std::size_t previous_num_particles = 0,
328 const uint64_t seed = 342343901 )
329{
330 using exec_space = typename PositionType::memory_space::execution_space;
331 createParticles( tag, exec_space{}, positions, num_particles, box_min,
332 box_max, previous_num_particles, seed );
333}
334
335//---------------------------------------------------------------------------//
336// Random with minimum separation
337//---------------------------------------------------------------------------//
338
367template <class ExecutionSpace, class InitFunctor, class ParticleListType,
368 class PositionTag,
369 template <class, std::size_t, class...> class ArrayType, class Scalar,
370 std::size_t NumSpaceDim, class... Args>
371int createParticles(
372 InitRandom tag, ExecutionSpace exec_space,
373 const InitFunctor& create_functor, ParticleListType& particle_list,
374 PositionTag position_tag, const std::size_t num_particles,
375 const double min_dist,
376#ifdef KOKKOS_ENABLE_DEPRECATED_CODE_4
377 const ArrayType<Scalar, NumSpaceDim, Args...> box_min,
378 const ArrayType<Scalar, NumSpaceDim, Args...> box_max,
379#else
380 const ArrayType<Scalar, NumSpaceDim> box_min,
381 const ArrayType<Scalar, NumSpaceDim> box_max,
382#endif
383 const std::size_t previous_num_particles = 0,
384 const bool shrink_to_fit = true, const uint64_t seed = 342343901,
385 typename std::enable_if<is_particle_list<ParticleListType>::value,
386 int>::type* = 0 )
387{
388 double min_dist_sqr = min_dist * min_dist;
389
390 // Resize the aosoa prior to lambda capture.
391 auto& aosoa = particle_list.aosoa();
392 aosoa.resize( previous_num_particles + num_particles );
393
394 auto positions = particle_list.slice( position_tag );
395
396 // Create the functor which ignores particles within the radius of another.
397 auto min_distance_op =
398 KOKKOS_LAMBDA( const int id, const double px[NumSpaceDim], const double,
399 typename ParticleListType::particle_type& particle )
400 {
401 // Ensure this particle is not within the minimum distance of any other
402 // existing particle.
403 for ( int n = 0; n < id; n++ )
404 {
405 double dist = 0.0;
406 for ( std::size_t d = 0; d < NumSpaceDim; ++d )
407 {
408 double dx = positions( n, d ) - px[d];
409 dist += dx * dx;
410 }
411 if ( dist < min_dist_sqr )
412 return false;
413 }
414
415 bool create = create_functor( id, px, 0.0, particle );
416 return create;
417 };
418
419 // Pass the functor to the general case.
420 return createParticles( tag, exec_space, min_distance_op, particle_list,
421 num_particles, box_min, box_max,
422 previous_num_particles, shrink_to_fit, seed );
423}
424
452template <class InitFunctor, class ParticleListType, class PositionTag,
453 template <class, std::size_t, class...> class ArrayType, class Scalar,
454 std::size_t NumSpaceDim, class... Args>
455int createParticles( InitRandom tag, const InitFunctor& create_functor,
456 ParticleListType& particle_list, PositionTag position_tag,
457 const std::size_t num_particles, const double min_dist,
458#ifdef KOKKOS_ENABLE_DEPRECATED_CODE_4
459 const ArrayType<Scalar, NumSpaceDim, Args...> box_min,
460 const ArrayType<Scalar, NumSpaceDim, Args...> box_max,
461#else
462 const ArrayType<Scalar, NumSpaceDim> box_min,
463 const ArrayType<Scalar, NumSpaceDim> box_max,
464#endif
465 const std::size_t previous_num_particles = 0,
466 const bool shrink_to_fit = true,
467 const uint64_t seed = 342343901 )
468{
469 using exec_space = typename ParticleListType::memory_space::execution_space;
470 return createParticles( tag, exec_space{}, create_functor, particle_list,
471 position_tag, num_particles, min_dist, box_min,
472 box_max, previous_num_particles, shrink_to_fit,
473 seed );
474}
475
491template <class ExecutionSpace, class PositionType, class Scalar>
492void createParticles( InitRandom tag, ExecutionSpace exec_space,
493 PositionType& positions, const std::size_t num_particles,
494 const double min_dist, const Scalar box_min[3],
495 const Scalar box_max[3],
496 const std::size_t previous_num_particles = 0,
497 const uint64_t seed = 342343901 )
498{
499 auto kokkos_min = copyArray<Scalar, 3>( box_min );
500 auto kokkos_max = copyArray<Scalar, 3>( box_max );
501 createParticles( tag, exec_space, positions, num_particles, min_dist,
502 kokkos_min, kokkos_max, previous_num_particles, seed );
503}
504
519template <class PositionType, class Scalar>
520void createParticles( InitRandom tag, PositionType& positions,
521 const std::size_t num_particles, const double min_dist,
522 const Scalar box_min[3], const Scalar box_max[3],
523 const std::size_t previous_num_particles = 0,
524 const uint64_t seed = 342343901 )
525{
526 using exec_space = typename PositionType::memory_space::execution_space;
527 createParticles( tag, exec_space{}, positions, num_particles, min_dist,
528 box_min, box_max, previous_num_particles, seed );
529}
530
531} // namespace Cabana
532
533#endif
Cabana_ParticleList.hpp
Application-level particle storage and single particle access.
Cabana_Slice.hpp
Slice a single particle property from an AoSoA.
Cabana_Utils.hpp
Cabana utilities.
Cabana
Core: particle data structures and algorithms.
Definition Cabana_AoSoA.hpp:36
Cabana::copyArray
auto copyArray(const std::array< Scalar, Dim > input)
Copy std::array into Kokkos::Array for potential device use.
Definition Cabana_Utils.hpp:32
Cabana::createParticles
int createParticles(InitRandom, ExecutionSpace exec_space, const InitFunctor &create_functor, ParticleListType &particle_list, const std::size_t num_particles, const ArrayType< Scalar, NumSpaceDim > box_min, const ArrayType< Scalar, NumSpaceDim > box_max, const std::size_t previous_num_particles=0, const bool shrink_to_fit=true, const uint64_t seed=342343901, typename std::enable_if< is_particle_list< ParticleListType >::value, int >::type *=0)
Initialize random particles given an initialization functor.
Definition Cabana_ParticleInit.hpp:74
Cabana::checkSize
void checkSize(SliceType slice, const std::size_t size, typename std::enable_if< is_slice< SliceType >::value, int >::type *=0)
Check slice size (differs from Kokkos View).
Definition Cabana_Slice.hpp:997
Cabana::InitRandom
Random particle initialization type tag.
Definition Cabana_ParticleInit.hpp:41
Cabana::InitUniform
Uniform particle initialization type tag.
Definition Cabana_ParticleInit.hpp:37
Cabana::is_particle_list
ParticleList static type checker.
Definition Cabana_ParticleList.hpp:282
Cabana::is_slice
Slice static type checker.
Definition Cabana_Slice.hpp:868