From 520bb2fc41a9b367eecca35d321a96ff5791a5e0 Mon Sep 17 00:00:00 2001 From: "Chong, Kwitae" Date: Mon, 23 Jan 2023 11:28:27 -0500 Subject: [PATCH 1/5] add unifor surface particle initialization --- src/Picasso_ParticleInit.hpp | 139 +++++++++++++++++++++++++++++++++++ 1 file changed, 139 insertions(+) diff --git a/src/Picasso_ParticleInit.hpp b/src/Picasso_ParticleInit.hpp index 30ef4d1..577b8de 100644 --- a/src/Picasso_ParticleInit.hpp +++ b/src/Picasso_ParticleInit.hpp @@ -535,6 +535,145 @@ void initializeParticlesSurface( InitRandom, const ExecutionSpace&, Kokkos::Profiling::popRegion(); } +//---------------------------------------------------------------------------// +/*! + \brief Initialize a uniform number of particles in each cell given an + initialization functor. + + \tparam ParticleListType The type of particle list to initialize. + + \tparam InitFunc Initialization functor type. See the documentation below + for the create_functor parameter on the signature of this functor. + + \tparam FacetGeometry The type of the facet geometry representing the + surface to initialize with particles + + \param Initialization type tag. + + \param particles_per_facet The number of particles to populate each cell + dimension with. + particles_per_facet = 1 : on the centroid + particles_per_facet = 3 : on the centers between centroid and vertex + particles_per_facet = 4 : on the centroid + + on the centers between centroid and vertex + + \param surface A FacetGeometry type (FacetGeoemtry or MarchingCubes data) + containing the list of surface facets to seed the particle list + + \param create_functor A functor which populates a particle on a facet 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], const double pan[3], double area, + typename ParticleAoSoA::tuple_type& particle ); + + px[3] : particle logical coordinates + pan[3] : facet area normal vector + area : the discretized surface area + particle& : the surface particle to create + + \param surface_particle_list The list of particles to populate. This will be + filled with particles and resized to a size equal to the number of particles + created. +*/ +template +void initializeParticlesSurface( InitUniform, const ExecutionSpace&, + const int particles_per_facet, + const FacetGeometry& surface, + const InitFunc& create_functor, + ParticleListType& surface_particle_list ) +{ + Kokkos::Profiling::pushRegion( "Picasso::initializeParticles::Surface" ); + + // Particle type. + using particle_type = typename ParticleListType::particle_type; + + // Get the local grid. + const auto& local_grid = *( surface_particle_list.mesh().localGrid() ); + + // Create a local mesh. + auto local_mesh = Cajita::createLocalMesh( local_grid ); + + // Get the global grid. + const auto& global_grid = local_grid.globalGrid(); + + // Get the particles. + auto& particles = surface_particle_list.aosoa(); + + // Allocate the particles. + int num_facets = surface.facets.extent( 0 ); + int num_particles = particles_per_facet * num_facets; + particles.resize( num_particles ); + + // FIXME: Re-thread over particles instead of facets. + // Initialize particles. + Kokkos::parallel_for( + "Picasso::ParticleInit::RandomSurface", + Kokkos::RangePolicy( 0, num_facets ), + KOKKOS_LAMBDA( const int f ) { + // Particle coordinate. + Vec3 px; + + // Particle. + particle_type particle; + + Vec3 a = { surface.facets( f, 0, 0 ), + surface.facets( f, 0, 1 ), + surface.facets( f, 0, 2 ) }; + + Vec3 b = { surface.facets( f, 1, 0 ), + surface.facets( f, 1, 1 ), + surface.facets( f, 1, 2 ) }; + + Vec3 c = { surface.facets( f, 2, 0 ), + surface.facets( f, 2, 1 ), + surface.facets( f, 2, 2 ) }; + + Vec3 ab = b - a; + Vec3 ac = c - a; + + auto cross = ab % ac; + double c2 = ~cross * cross; + double sqc2 = sqrt( c2 ); + double facet_area = 0.5 * sqc2; + + Vec3 pan = + cross / sqc2; // Create unit normal vector for the facet area + + // Particle surface area. + double pa = facet_area / particles_per_facet; + + // centroid of triangle facet + Vec3 centroid = ( a + b + c ) / 3.0; + Vec3 a_centroid = ( a + centroid ) / 2.0; + Vec3 b_centroid = ( b + centroid ) / 2.0; + Vec3 c_centroid = ( c + centroid ) / 2.0; + Kokkos::Array, 4> abc_centroid + = { centroid, a_centroid, b_centroid, c_centroid }; + + for ( int ip = 0; ip < particles_per_facet; ++ip ) + { + // Local particle id. + int pid = f * particles_per_facet + ip; + + // ppf = 3 + if( particles_per_facet == 3 ) + px = abc_centroid[ip+1]; + // ppf = 1 or 4 + else + px = abc_centroid[ip]; + + // Create a new particle with the given logical + // coordinates. + create_functor( px.data(), pan.data(), pa, particle ); + + particles.setTuple( pid, particle.tuple() ); + } + } ); + + Kokkos::Profiling::popRegion(); +} //---------------------------------------------------------------------------// } // end namespace Picasso From e6981bc15877cb5bc59d83a97f581e7b4815d1fb Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Thu, 26 Jan 2023 17:21:09 -0500 Subject: [PATCH 2/5] format --- src/Picasso_ParticleInit.hpp | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/Picasso_ParticleInit.hpp b/src/Picasso_ParticleInit.hpp index 577b8de..78e7931 100644 --- a/src/Picasso_ParticleInit.hpp +++ b/src/Picasso_ParticleInit.hpp @@ -646,11 +646,11 @@ void initializeParticlesSurface( InitUniform, const ExecutionSpace&, // centroid of triangle facet Vec3 centroid = ( a + b + c ) / 3.0; - Vec3 a_centroid = ( a + centroid ) / 2.0; - Vec3 b_centroid = ( b + centroid ) / 2.0; - Vec3 c_centroid = ( c + centroid ) / 2.0; - Kokkos::Array, 4> abc_centroid - = { centroid, a_centroid, b_centroid, c_centroid }; + Vec3 a_centroid = ( a + centroid ) / 2.0; + Vec3 b_centroid = ( b + centroid ) / 2.0; + Vec3 c_centroid = ( c + centroid ) / 2.0; + Kokkos::Array, 4> abc_centroid = { + centroid, a_centroid, b_centroid, c_centroid }; for ( int ip = 0; ip < particles_per_facet; ++ip ) { @@ -658,11 +658,11 @@ void initializeParticlesSurface( InitUniform, const ExecutionSpace&, int pid = f * particles_per_facet + ip; // ppf = 3 - if( particles_per_facet == 3 ) - px = abc_centroid[ip+1]; - // ppf = 1 or 4 - else - px = abc_centroid[ip]; + if ( particles_per_facet == 3 ) + px = abc_centroid[ip + 1]; + // ppf = 1 or 4 + else + px = abc_centroid[ip]; // Create a new particle with the given logical // coordinates. From ea74a236b6cd31fa33cf2ba21fbe16df5f847f1c Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 27 Jan 2023 12:39:37 -0500 Subject: [PATCH 3/5] Input surface particles per facet median to match internal uniform style --- src/Picasso_ParticleInit.hpp | 62 +++++++++++++++++------------------- 1 file changed, 30 insertions(+), 32 deletions(-) diff --git a/src/Picasso_ParticleInit.hpp b/src/Picasso_ParticleInit.hpp index 78e7931..dba1f88 100644 --- a/src/Picasso_ParticleInit.hpp +++ b/src/Picasso_ParticleInit.hpp @@ -550,12 +550,8 @@ void initializeParticlesSurface( InitRandom, const ExecutionSpace&, \param Initialization type tag. - \param particles_per_facet The number of particles to populate each cell - dimension with. - particles_per_facet = 1 : on the centroid - particles_per_facet = 3 : on the centers between centroid and vertex - particles_per_facet = 4 : on the centroid + - on the centers between centroid and vertex + \param particles_per_facet_median The number of particles to populate each + facet median with (+1 on the centroid for 3n+1 total per facet). \param surface A FacetGeometry type (FacetGeoemtry or MarchingCubes data) containing the list of surface facets to seed the particle list @@ -579,7 +575,7 @@ void initializeParticlesSurface( InitRandom, const ExecutionSpace&, template void initializeParticlesSurface( InitUniform, const ExecutionSpace&, - const int particles_per_facet, + const int particles_per_facet_median, const FacetGeometry& surface, const InitFunc& create_functor, ParticleListType& surface_particle_list ) @@ -603,6 +599,7 @@ void initializeParticlesSurface( InitUniform, const ExecutionSpace&, // Allocate the particles. int num_facets = surface.facets.extent( 0 ); + int particles_per_facet = particles_per_facet_median * 3 + 1; int num_particles = particles_per_facet * num_facets; particles.resize( num_particles ); @@ -638,38 +635,39 @@ void initializeParticlesSurface( InitUniform, const ExecutionSpace&, double sqc2 = sqrt( c2 ); double facet_area = 0.5 * sqc2; - Vec3 pan = - cross / sqc2; // Create unit normal vector for the facet area - + // Create unit normal vector for the facet area + Vec3 pan = cross / sqc2; // Particle surface area. double pa = facet_area / particles_per_facet; // centroid of triangle facet Vec3 centroid = ( a + b + c ) / 3.0; - Vec3 a_centroid = ( a + centroid ) / 2.0; - Vec3 b_centroid = ( b + centroid ) / 2.0; - Vec3 c_centroid = ( c + centroid ) / 2.0; - Kokkos::Array, 4> abc_centroid = { - centroid, a_centroid, b_centroid, c_centroid }; - - for ( int ip = 0; ip < particles_per_facet; ++ip ) - { - // Local particle id. - int pid = f * particles_per_facet + ip; + int median_div = particles_per_facet_median + 2; + Vec3 a_median = ( a + centroid ) / median_div; + Vec3 b_median = ( b + centroid ) / median_div; + Vec3 c_median = ( c + centroid ) / median_div; + Kokkos::Array, 3> abc_median = { a_median, b_median, + c_median }; + + // Create one particle at the centroid. + create_functor( centroid.data(), pan.data(), pa, particle ); + particles.setTuple( f * particles_per_facet, particle.tuple() ); + + int ppf_count = 1; + for ( int it = 0; it < 3; ++it ) + for ( int ip = 0; ip < particles_per_facet_median; ++ip ) + { + // Local particle id. + int pid = f * particles_per_facet + ppf_count; - // ppf = 3 - if ( particles_per_facet == 3 ) - px = abc_centroid[ip + 1]; - // ppf = 1 or 4 - else - px = abc_centroid[ip]; + px = abc_median[it]; + // abc_median[it] = ip * abc_median[it]; - // Create a new particle with the given logical - // coordinates. - create_functor( px.data(), pan.data(), pa, particle ); - - particles.setTuple( pid, particle.tuple() ); - } + // Create a new particle with the given logical coordinates. + create_functor( px.data(), pan.data(), pa, particle ); + particles.setTuple( pid, particle.tuple() ); + ppf_count++; + } } ); Kokkos::Profiling::popRegion(); From e33a94eb1fe23ff5588308b8ec51b08e61c4a0c4 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 27 Jan 2023 12:41:40 -0500 Subject: [PATCH 4/5] Update labels --- src/Picasso_ParticleInit.hpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/Picasso_ParticleInit.hpp b/src/Picasso_ParticleInit.hpp index dba1f88..d53571d 100644 --- a/src/Picasso_ParticleInit.hpp +++ b/src/Picasso_ParticleInit.hpp @@ -438,7 +438,8 @@ void initializeParticlesSurface( InitRandom, const ExecutionSpace&, const InitFunc& create_functor, ParticleListType& surface_particle_list ) { - Kokkos::Profiling::pushRegion( "Picasso::initializeParticles::Surface" ); + Kokkos::Profiling::pushRegion( + "Picasso::initializeParticles::RandomSurface" ); // Particle type. using particle_type = typename ParticleListType::particle_type; @@ -580,7 +581,8 @@ void initializeParticlesSurface( InitUniform, const ExecutionSpace&, const InitFunc& create_functor, ParticleListType& surface_particle_list ) { - Kokkos::Profiling::pushRegion( "Picasso::initializeParticles::Surface" ); + Kokkos::Profiling::pushRegion( + "Picasso::initializeParticles::UniformSurface" ); // Particle type. using particle_type = typename ParticleListType::particle_type; @@ -606,7 +608,7 @@ void initializeParticlesSurface( InitUniform, const ExecutionSpace&, // FIXME: Re-thread over particles instead of facets. // Initialize particles. Kokkos::parallel_for( - "Picasso::ParticleInit::RandomSurface", + "Picasso::ParticleInit::UniformSurface", Kokkos::RangePolicy( 0, num_facets ), KOKKOS_LAMBDA( const int f ) { // Particle coordinate. From 58de92a6b00433f242ec070314c539e498c9943f Mon Sep 17 00:00:00 2001 From: "Chong, Kwitae" Date: Sat, 28 Jan 2023 13:02:45 -0500 Subject: [PATCH 5/5] fix uniform surface particle to be generated on the internal division points by particles_per_facet_median --- src/Picasso_ParticleInit.hpp | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/Picasso_ParticleInit.hpp b/src/Picasso_ParticleInit.hpp index d53571d..bef3a5b 100644 --- a/src/Picasso_ParticleInit.hpp +++ b/src/Picasso_ParticleInit.hpp @@ -644,21 +644,24 @@ void initializeParticlesSurface( InitUniform, const ExecutionSpace&, // centroid of triangle facet Vec3 centroid = ( a + b + c ) / 3.0; - int median_div = particles_per_facet_median + 2; - Vec3 a_median = ( a + centroid ) / median_div; - Vec3 b_median = ( b + centroid ) / median_div; - Vec3 c_median = ( c + centroid ) / median_div; - Kokkos::Array, 3> abc_median = { a_median, b_median, - c_median }; + int median_div = particles_per_facet_median + 1; // Create one particle at the centroid. create_functor( centroid.data(), pan.data(), pa, particle ); particles.setTuple( f * particles_per_facet, particle.tuple() ); int ppf_count = 1; - for ( int it = 0; it < 3; ++it ) - for ( int ip = 0; ip < particles_per_facet_median; ++ip ) + for ( int ip = 1; ip <= particles_per_facet_median; ++ip ) + { + for ( int it = 0; it < 3; ++it ) { + // internal division points + Vec3 a_median = ( a*(median_div-ip) + centroid*ip ) / median_div; + Vec3 b_median = ( b*(median_div-ip) + centroid*ip ) / median_div; + Vec3 c_median = ( c*(median_div-ip) + centroid*ip ) / median_div; + Kokkos::Array, 3> abc_median = { a_median, b_median, + c_median }; + // Local particle id. int pid = f * particles_per_facet + ppf_count; @@ -670,6 +673,7 @@ void initializeParticlesSurface( InitUniform, const ExecutionSpace&, particles.setTuple( pid, particle.tuple() ); ppf_count++; } + } } ); Kokkos::Profiling::popRegion();