From b50761ce2a41a5362428d8b0aa9a055da9288b1e Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 3 Jul 2024 15:35:17 -0400 Subject: [PATCH 1/6] Rewrite boundary conditions for multiple geometries --- examples/mechanics/crack_branching.cpp | 15 ++-- examples/mechanics/kalthoff_winkler.cpp | 5 +- src/CabanaPD_Boundary.hpp | 91 +++++++++++++++++++------ 3 files changed, 80 insertions(+), 31 deletions(-) diff --git a/examples/mechanics/crack_branching.cpp b/examples/mechanics/crack_branching.cpp index b9dc1b7d..6833c5fd 100644 --- a/examples/mechanics/crack_branching.cpp +++ b/examples/mechanics/crack_branching.cpp @@ -90,13 +90,14 @@ void crackBranchingExample( const std::string filename ) double dy = particles->dx[1]; double b0 = sigma0 / dy; - CabanaPD::RegionBoundary plane1( low_corner[0], high_corner[0], - low_corner[1] - dy, low_corner[1] + dy, - low_corner[2], high_corner[2] ); - CabanaPD::RegionBoundary plane2( low_corner[0], high_corner[0], - high_corner[1] - dy, high_corner[1] + dy, - low_corner[2], high_corner[2] ); - std::vector planes = { plane1, plane2 }; + CabanaPD::RegionBoundary plane1( + low_corner[0], high_corner[0], low_corner[1] - dy, low_corner[1] + dy, + low_corner[2], high_corner[2] ); + CabanaPD::RegionBoundary plane2( + low_corner[0], high_corner[0], high_corner[1] - dy, high_corner[1] + dy, + low_corner[2], high_corner[2] ); + std::vector> planes = { + plane1, plane2 }; auto particles_f = particles->getForce(); auto particles_x = particles->getReferencePosition(); // Create a symmetric force BC in the y-direction. diff --git a/examples/mechanics/kalthoff_winkler.cpp b/examples/mechanics/kalthoff_winkler.cpp index 9da84ce1..69a3a133 100644 --- a/examples/mechanics/kalthoff_winkler.cpp +++ b/examples/mechanics/kalthoff_winkler.cpp @@ -98,10 +98,11 @@ void kalthoffWinklerExample( const std::string filename ) // ==================================================== double dx = particles->dx[0]; double x_bc = -0.5 * height; - CabanaPD::RegionBoundary plane( + CabanaPD::RegionBoundary plane( x_bc - dx, x_bc + dx * 1.25, y_prenotch1 - dx * 0.25, y_prenotch2 + dx * 0.25, -thickness, thickness ); - std::vector planes = { plane }; + std::vector> planes = { + plane }; auto bc = createBoundaryCondition( CabanaPD::ForceValueBCTag{}, 0.0, exec_space{}, *particles, planes ); diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp index 5ac27da0..3eedd7b2 100644 --- a/src/CabanaPD_Boundary.hpp +++ b/src/CabanaPD_Boundary.hpp @@ -21,8 +21,20 @@ namespace CabanaPD { -// Define a plane or other rectilinear subset of the system as the boundary. -struct RegionBoundary +struct RectangularPrism +{ +}; +struct Cylinder +{ +}; + +// Forward declaration. +template +struct RegionBoundary; + +// Define a subset of the system as the boundary with a rectangular prism. +template <> +struct RegionBoundary { double low_x; double high_x; @@ -40,14 +52,53 @@ struct RegionBoundary , high_y( _high_y ) , low_z( _low_z ) , high_z( _high_z ){}; + + template + KOKKOS_INLINE_FUNCTION bool inside( const PositionType& x, + const int pid ) const + { + return ( x( pid, 0 ) >= low_x && x( pid, 0 ) <= high_x && + x( pid, 1 ) >= low_y && x( pid, 1 ) <= high_y && + x( pid, 2 ) >= low_z && x( pid, 2 ) <= high_z ); + } +}; + +// Define a subset of the system as the boundary with a cylinder. +template <> +struct RegionBoundary +{ + double radius_in; + double radius_out; + double low_z; + double high_z; + double x_center = 0.0; + double y_center = 0.0; + + RegionBoundary( const double _radius_in, const double _radius_out, + const double _low_z, const double _high_z ) + : radius_in( _radius_in ) + , radius_out( _radius_out ) + , low_z( _low_z ) + , high_z( _high_z ){}; + + template + KOKKOS_INLINE_FUNCTION bool inside( const PositionType& x, + const int pid ) const + { + double rsq = ( x( pid, 0 ) - x_center ) * ( x( pid, 0 ) - x_center ) + + ( x( pid, 1 ) - y_center ) * ( x( pid, 1 ) - y_center ); + return ( rsq >= radius_in * radius_in && + rsq <= radius_out * radius_out && x( pid, 2 ) >= low_z && + x( pid, 2 ) <= high_z ); + } }; template struct BoundaryIndexSpace; // FIXME: fails for some cases if initial guess is not sufficient. -template -struct BoundaryIndexSpace +template +struct BoundaryIndexSpace> { using index_view_type = Kokkos::View; index_view_type _view; @@ -60,7 +111,7 @@ struct BoundaryIndexSpace template BoundaryIndexSpace( ExecSpace exec_space, Particles particles, - std::vector planes, + std::vector> planes, const double initial_guess ) { _timer.start(); @@ -85,7 +136,8 @@ struct BoundaryIndexSpace } template - void update( ExecSpace, Particles particles, RegionBoundary plane ) + void update( ExecSpace, Particles particles, + RegionBoundary region ) { auto count_host = Kokkos::create_mirror_view_and_copy( Kokkos::HostSpace{}, _count ); @@ -97,9 +149,7 @@ struct BoundaryIndexSpace Kokkos::RangePolicy policy( 0, particles.n_local ); auto index_functor = KOKKOS_LAMBDA( const std::size_t pid ) { - if ( x( pid, 0 ) >= plane.low_x && x( pid, 0 ) <= plane.high_x && - x( pid, 1 ) >= plane.low_y && x( pid, 1 ) <= plane.high_y && - x( pid, 2 ) >= plane.low_z && x( pid, 2 ) <= plane.high_z ) + if ( region.inside( x, pid ) ) { // Resize after count if needed. auto c = Kokkos::atomic_fetch_add( &count( 0 ), 1 ); @@ -124,9 +174,9 @@ struct BoundaryIndexSpace auto time() { return _timer.time(); }; }; -template +template auto createBoundaryIndexSpace( ExecSpace exec_space, Particles particles, - std::vector planes, + std::vector planes, const double initial_guess ) { using memory_space = typename Particles::memory_space; @@ -159,9 +209,8 @@ struct BoundaryCondition { } - template - void update( ExecSpace exec_space, Particles particles, - RegionBoundary plane ) + template + void update( ExecSpace exec_space, Particles particles, BoundaryType plane ) { _index_space.update( exec_space, particles, plane ); } @@ -202,9 +251,8 @@ struct BoundaryCondition { } - template - void update( ExecSpace exec_space, Particles particles, - RegionBoundary plane ) + template + void update( ExecSpace exec_space, Particles particles, BoundaryType plane ) { _index_space.update( exec_space, particles, plane ); } @@ -247,9 +295,8 @@ struct BoundaryCondition { } - template - void update( ExecSpace exec_space, Particles particles, - RegionBoundary plane ) + template + void update( ExecSpace exec_space, Particles particles, BoundaryType plane ) { _index_space.update( exec_space, particles, plane ); } @@ -288,7 +335,7 @@ auto createBoundaryCondition( BCTag, const double value, ExecSpace exec_space, { using memory_space = typename Particles::memory_space; using bc_index_type = BoundaryIndexSpace; - bc_index_type bc_indices = createBoundaryIndexSpace( + bc_index_type bc_indices = createBoundaryIndexSpace( exec_space, particles, planes, initial_guess ); return BoundaryCondition( value, bc_indices ); } @@ -304,7 +351,7 @@ auto createBoundaryCondition( UserFunctor user_functor, ExecSpace exec_space, { using memory_space = typename Particles::memory_space; using bc_index_type = BoundaryIndexSpace; - bc_index_type bc_indices = createBoundaryIndexSpace( + bc_index_type bc_indices = createBoundaryIndexSpace( exec_space, particles, planes, initial_guess ); return BoundaryCondition( bc_indices, user_functor, force_update ); From 111f2faefe93b4679313d95f1c545bcf33cb38b3 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 3 Jul 2024 15:39:02 -0400 Subject: [PATCH 2/6] Add BC convenience interface for single plane --- examples/mechanics/kalthoff_winkler.cpp | 4 +--- src/CabanaPD_Boundary.hpp | 23 +++++++++++++++++++++++ 2 files changed, 24 insertions(+), 3 deletions(-) diff --git a/examples/mechanics/kalthoff_winkler.cpp b/examples/mechanics/kalthoff_winkler.cpp index 69a3a133..a197c279 100644 --- a/examples/mechanics/kalthoff_winkler.cpp +++ b/examples/mechanics/kalthoff_winkler.cpp @@ -101,10 +101,8 @@ void kalthoffWinklerExample( const std::string filename ) CabanaPD::RegionBoundary plane( x_bc - dx, x_bc + dx * 1.25, y_prenotch1 - dx * 0.25, y_prenotch2 + dx * 0.25, -thickness, thickness ); - std::vector> planes = { - plane }; auto bc = createBoundaryCondition( CabanaPD::ForceValueBCTag{}, 0.0, - exec_space{}, *particles, planes ); + exec_space{}, *particles, plane ); // ==================================================== // Custom particle initialization diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp index 3eedd7b2..6d2cffcb 100644 --- a/src/CabanaPD_Boundary.hpp +++ b/src/CabanaPD_Boundary.hpp @@ -357,6 +357,29 @@ auto createBoundaryCondition( UserFunctor user_functor, ExecSpace exec_space, bc_indices, user_functor, force_update ); } +template +auto createBoundaryCondition( BCTag tag, const double value, + ExecSpace exec_space, Particles particles, + BoundaryType plane, + const double initial_guess = 0.5 ) +{ + std::vector plane_vec = { plane }; + return createBoundaryCondition( tag, value, exec_space, particles, + plane_vec, initial_guess ); +} + +template +auto createBoundaryCondition( UserFunctor user_functor, ExecSpace exec_space, + Particles particles, BoundaryType plane, + const bool force_update, + const double initial_guess = 0.5 ) +{ + std::vector plane_vec = { plane }; + return createBoundaryCondition( user_functor, exec_space, particles, + plane_vec, force_update, initial_guess ); +} + } // namespace CabanaPD #endif From 78ce23d560214b42af54933be6172744c0b46967 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Thu, 1 Aug 2024 11:49:31 -0400 Subject: [PATCH 3/6] Pass time in all BC (even if unused) --- examples/mechanics/crack_branching.cpp | 2 +- src/CabanaPD_Boundary.hpp | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/mechanics/crack_branching.cpp b/examples/mechanics/crack_branching.cpp index 6833c5fd..d34b2a50 100644 --- a/examples/mechanics/crack_branching.cpp +++ b/examples/mechanics/crack_branching.cpp @@ -101,7 +101,7 @@ void crackBranchingExample( const std::string filename ) auto particles_f = particles->getForce(); auto particles_x = particles->getReferencePosition(); // Create a symmetric force BC in the y-direction. - auto bc_op = KOKKOS_LAMBDA( const int pid ) + auto bc_op = KOKKOS_LAMBDA( const int pid, const double ) { // Get a modifiable copy of force. auto p_f = particles_f.getParticleView( pid ); diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp index 6d2cffcb..715144b2 100644 --- a/src/CabanaPD_Boundary.hpp +++ b/src/CabanaPD_Boundary.hpp @@ -216,7 +216,7 @@ struct BoundaryCondition } template - void apply( ExecSpace, ParticleType&, double ) + void apply( ExecSpace, ParticleType&, double time ) { _timer.start(); auto user = _user_functor; @@ -225,7 +225,7 @@ struct BoundaryCondition Kokkos::parallel_for( "CabanaPD::BC::apply", policy, KOKKOS_LAMBDA( const int b ) { auto pid = index_space( b ); - user( pid ); + user( pid, time ); } ); _timer.stop(); } From 3d5e2bb689d650d0c1ff66e89951ec65a196e0e3 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 7 Aug 2024 17:42:16 -0400 Subject: [PATCH 4/6] Add option for custom boundary index spaces --- src/CabanaPD_Boundary.hpp | 48 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp index 715144b2..26531b37 100644 --- a/src/CabanaPD_Boundary.hpp +++ b/src/CabanaPD_Boundary.hpp @@ -27,6 +27,9 @@ struct RectangularPrism struct Cylinder { }; +struct Custom +{ +}; // Forward declaration. template @@ -109,6 +112,7 @@ struct BoundaryIndexSpace> // Default for empty case. BoundaryIndexSpace() {} + // Construct from region (search for boundary particles). template BoundaryIndexSpace( ExecSpace exec_space, Particles particles, std::vector> planes, @@ -174,6 +178,19 @@ struct BoundaryIndexSpace> auto time() { return _timer.time(); }; }; +template +struct BoundaryIndexSpace +{ + using index_view_type = Kokkos::View; + index_view_type _view; + + // Construct from a View of boundary particles. + BoundaryIndexSpace( index_view_type input_view ) + : _view( input_view ) + { + } +}; + template auto createBoundaryIndexSpace( ExecSpace exec_space, Particles particles, std::vector planes, @@ -184,6 +201,13 @@ auto createBoundaryIndexSpace( ExecSpace exec_space, Particles particles, exec_space, particles, planes, initial_guess ); } +template +auto createBoundaryIndexSpace( BoundaryParticles particles ) +{ + using memory_space = typename BoundaryParticles::memory_space; + return BoundaryIndexSpace( particles ); +} + struct ForceValueBCTag { }; @@ -357,6 +381,30 @@ auto createBoundaryCondition( UserFunctor user_functor, ExecSpace exec_space, bc_indices, user_functor, force_update ); } +// Custom index space cases +template +auto createBoundaryCondition( BCTag, const double value, + BoundaryParticles particles ) +{ + using memory_space = typename BoundaryParticles::memory_space; + using bc_index_type = BoundaryIndexSpace; + bc_index_type bc_indices = createBoundaryIndexSpace( particles ); + return BoundaryCondition( value, bc_indices ); +} + +template +auto createBoundaryCondition( UserFunctor user_functor, + BoundaryParticles particles, + const bool force_update ) +{ + using memory_space = typename BoundaryParticles::memory_space; + using bc_index_type = BoundaryIndexSpace; + bc_index_type bc_indices = createBoundaryIndexSpace( particles ); + return BoundaryCondition( + bc_indices, user_functor, force_update ); +} + +// Wrappers for single plane cases. template auto createBoundaryCondition( BCTag tag, const double value, ExecSpace exec_space, Particles particles, From 92e877657d5327ddb337c9a177be02d31435d94d Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Thu, 15 Aug 2024 18:01:23 -0400 Subject: [PATCH 5/6] fixup: add center to the region constructor --- src/CabanaPD_Boundary.hpp | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp index 26531b37..574a8f95 100644 --- a/src/CabanaPD_Boundary.hpp +++ b/src/CabanaPD_Boundary.hpp @@ -74,15 +74,18 @@ struct RegionBoundary double radius_out; double low_z; double high_z; - double x_center = 0.0; - double y_center = 0.0; + double x_center; + double y_center; RegionBoundary( const double _radius_in, const double _radius_out, - const double _low_z, const double _high_z ) + const double _low_z, const double _high_z, + double _x_center = 0.0, double _y_center = 0.0 ) : radius_in( _radius_in ) , radius_out( _radius_out ) , low_z( _low_z ) - , high_z( _high_z ){}; + , high_z( _high_z ) + , x_center( _x_center ) + , y_center( _y_center ){}; template KOKKOS_INLINE_FUNCTION bool inside( const PositionType& x, From 95908c335a9d19d871dbde930501f62068c0009b Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 16 Aug 2024 10:46:07 -0400 Subject: [PATCH 6/6] fixup: add custom region BC option --- src/CabanaPD_Boundary.hpp | 22 +++++++++++++++++++--- 1 file changed, 19 insertions(+), 3 deletions(-) diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp index 574a8f95..c1a3057c 100644 --- a/src/CabanaPD_Boundary.hpp +++ b/src/CabanaPD_Boundary.hpp @@ -31,9 +31,25 @@ struct Custom { }; -// Forward declaration. -template -struct RegionBoundary; +// User-specifed custom boundary. Must use the signature: +// bool operator()(PositionType&, const int) +template +struct RegionBoundary +{ + UserFunctor _user_functor; + + RegionBoundary( UserFunctor user ) + : _user_functor( user ) + { + } + + template + KOKKOS_INLINE_FUNCTION bool inside( const PositionType& x, + const int pid ) const + { + return user( x, pid ); + } +}; // Define a subset of the system as the boundary with a rectangular prism. template <>