Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Generalize boundary condition options #125

Merged
merged 6 commits into from
Aug 16, 2024
Merged

Generalize boundary condition options #125

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
b50761c
Rewrite boundary conditions for multiple geometries
streeve Jul 3, 2024
111f2fa
Add BC convenience interface for single plane
streeve Jul 3, 2024
78ce23d
Pass time in all BC (even if unused)
streeve Aug 1, 2024
3d5e2bb
Add option for custom boundary index spaces
streeve Aug 7, 2024
92e8776
fixup: add center to the region constructor
streeve Aug 15, 2024
95908c3
fixup: add custom region BC option
streeve Aug 16, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 9 additions & 8 deletions examples/mechanics/crack_branching.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -90,17 +90,18 @@ 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<CabanaPD::RegionBoundary> planes = { plane1, plane2 };
CabanaPD::RegionBoundary<CabanaPD::RectangularPrism> plane1(
low_corner[0], high_corner[0], low_corner[1] - dy, low_corner[1] + dy,
low_corner[2], high_corner[2] );
CabanaPD::RegionBoundary<CabanaPD::RectangularPrism> plane2(
low_corner[0], high_corner[0], high_corner[1] - dy, high_corner[1] + dy,
low_corner[2], high_corner[2] );
std::vector<CabanaPD::RegionBoundary<CabanaPD::RectangularPrism>> planes = {
streeve marked this conversation as resolved.
Show resolved Hide resolved
plane1, plane2 };
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 )
streeve marked this conversation as resolved.
Show resolved Hide resolved
{
// Get a modifiable copy of force.
auto p_f = particles_f.getParticleView( pid );
Expand Down
5 changes: 2 additions & 3 deletions examples/mechanics/kalthoff_winkler.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -98,12 +98,11 @@ void kalthoffWinklerExample( const std::string filename )
// ====================================================
double dx = particles->dx[0];
double x_bc = -0.5 * height;
CabanaPD::RegionBoundary plane(
CabanaPD::RegionBoundary<CabanaPD::RectangularPrism> plane(
x_bc - dx, x_bc + dx * 1.25, y_prenotch1 - dx * 0.25,
streeve marked this conversation as resolved.
Show resolved Hide resolved
y_prenotch2 + dx * 0.25, -thickness, thickness );
std::vector<CabanaPD::RegionBoundary> planes = { plane };
auto bc = createBoundaryCondition( CabanaPD::ForceValueBCTag{}, 0.0,
exec_space{}, *particles, planes );
exec_space{}, *particles, plane );

// ====================================================
// Custom particle initialization
Expand Down
183 changes: 160 additions & 23 deletions src/CabanaPD_Boundary.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,8 +21,39 @@
namespace CabanaPD
{

// Define a plane or other rectilinear subset of the system as the boundary.
struct RectangularPrism
{
};
struct Cylinder
{
};
struct Custom
{
};

// User-specifed custom boundary. Must use the signature:
// bool operator()(PositionType&, const int)
template <typename UserFunctor>
struct RegionBoundary
{
UserFunctor _user_functor;

RegionBoundary( UserFunctor user )
: _user_functor( user )
{
}

template <class PositionType>
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 <>
struct RegionBoundary<RectangularPrism>
{
double low_x;
double high_x;
Expand All @@ -40,14 +71,56 @@ struct RegionBoundary
, high_y( _high_y )
, low_z( _low_z )
, high_z( _high_z ){};

template <class PositionType>
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<Cylinder>
{
double radius_in;
double radius_out;
double low_z;
double high_z;
double x_center;
double y_center;

RegionBoundary( const double _radius_in, const double _radius_out,
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 )
, x_center( _x_center )
, y_center( _y_center ){};

template <class PositionType>
KOKKOS_INLINE_FUNCTION bool inside( const PositionType& x,
const int pid ) const
{
double rsq = ( x( pid, 0 ) - x_center ) * ( x( pid, 0 ) - x_center ) +
streeve marked this conversation as resolved.
Show resolved Hide resolved
( 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 <class MemorySpace, class BoundaryType>
struct BoundaryIndexSpace;

// FIXME: fails for some cases if initial guess is not sufficient.
template <class MemorySpace>
struct BoundaryIndexSpace<MemorySpace, RegionBoundary>
template <class MemorySpace, class GeometryType>
struct BoundaryIndexSpace<MemorySpace, RegionBoundary<GeometryType>>
{
using index_view_type = Kokkos::View<std::size_t*, MemorySpace>;
index_view_type _view;
Expand All @@ -58,9 +131,10 @@ struct BoundaryIndexSpace<MemorySpace, RegionBoundary>
// Default for empty case.
BoundaryIndexSpace() {}

// Construct from region (search for boundary particles).
template <class ExecSpace, class Particles>
BoundaryIndexSpace( ExecSpace exec_space, Particles particles,
std::vector<RegionBoundary> planes,
std::vector<RegionBoundary<GeometryType>> planes,
const double initial_guess )
{
_timer.start();
Expand All @@ -85,7 +159,8 @@ struct BoundaryIndexSpace<MemorySpace, RegionBoundary>
}

template <class ExecSpace, class Particles>
void update( ExecSpace, Particles particles, RegionBoundary plane )
void update( ExecSpace, Particles particles,
RegionBoundary<GeometryType> region )
{
auto count_host =
Kokkos::create_mirror_view_and_copy( Kokkos::HostSpace{}, _count );
Expand All @@ -97,9 +172,7 @@ struct BoundaryIndexSpace<MemorySpace, RegionBoundary>
Kokkos::RangePolicy<ExecSpace> 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 );
Expand All @@ -124,16 +197,36 @@ struct BoundaryIndexSpace<MemorySpace, RegionBoundary>
auto time() { return _timer.time(); };
};

template <class BoundaryType, class ExecSpace, class Particles>
template <class MemorySpace>
struct BoundaryIndexSpace<MemorySpace, Custom>
{
using index_view_type = Kokkos::View<std::size_t*, MemorySpace>;
index_view_type _view;

// Construct from a View of boundary particles.
BoundaryIndexSpace( index_view_type input_view )
: _view( input_view )
{
}
};

template <class ExecSpace, class Particles, class BoundaryType>
auto createBoundaryIndexSpace( ExecSpace exec_space, Particles particles,
std::vector<RegionBoundary> planes,
std::vector<BoundaryType> planes,
const double initial_guess )
{
using memory_space = typename Particles::memory_space;
return BoundaryIndexSpace<memory_space, BoundaryType>(
exec_space, particles, planes, initial_guess );
}

template <class BoundaryParticles>
auto createBoundaryIndexSpace( BoundaryParticles particles )
{
using memory_space = typename BoundaryParticles::memory_space;
return BoundaryIndexSpace<memory_space, Custom>( particles );
}

struct ForceValueBCTag
{
};
Expand All @@ -159,15 +252,14 @@ struct BoundaryCondition
{
}

template <class ExecSpace, class Particles>
void update( ExecSpace exec_space, Particles particles,
RegionBoundary plane )
template <class ExecSpace, class Particles, class BoundaryType>
void update( ExecSpace exec_space, Particles particles, BoundaryType plane )
{
_index_space.update( exec_space, particles, plane );
}

template <class ExecSpace, class ParticleType>
void apply( ExecSpace, ParticleType&, double )
void apply( ExecSpace, ParticleType&, double time )
streeve marked this conversation as resolved.
Show resolved Hide resolved
{
_timer.start();
auto user = _user_functor;
Expand All @@ -176,7 +268,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 );
streeve marked this conversation as resolved.
Show resolved Hide resolved
} );
_timer.stop();
}
Expand All @@ -202,9 +294,8 @@ struct BoundaryCondition<BCIndexSpace, ForceValueBCTag>
{
}

template <class ExecSpace, class Particles>
void update( ExecSpace exec_space, Particles particles,
RegionBoundary plane )
template <class ExecSpace, class Particles, class BoundaryType>
void update( ExecSpace exec_space, Particles particles, BoundaryType plane )
{
_index_space.update( exec_space, particles, plane );
}
Expand Down Expand Up @@ -247,9 +338,8 @@ struct BoundaryCondition<BCIndexSpace, ForceUpdateBCTag>
{
}

template <class ExecSpace, class Particles>
void update( ExecSpace exec_space, Particles particles,
RegionBoundary plane )
template <class ExecSpace, class Particles, class BoundaryType>
void update( ExecSpace exec_space, Particles particles, BoundaryType plane )
{
_index_space.update( exec_space, particles, plane );
}
Expand Down Expand Up @@ -288,7 +378,7 @@ auto createBoundaryCondition( BCTag, const double value, ExecSpace exec_space,
{
using memory_space = typename Particles::memory_space;
using bc_index_type = BoundaryIndexSpace<memory_space, BoundaryType>;
bc_index_type bc_indices = createBoundaryIndexSpace<BoundaryType>(
bc_index_type bc_indices = createBoundaryIndexSpace(
exec_space, particles, planes, initial_guess );
return BoundaryCondition<bc_index_type, BCTag>( value, bc_indices );
}
Expand All @@ -304,12 +394,59 @@ auto createBoundaryCondition( UserFunctor user_functor, ExecSpace exec_space,
{
using memory_space = typename Particles::memory_space;
using bc_index_type = BoundaryIndexSpace<memory_space, BoundaryType>;
bc_index_type bc_indices = createBoundaryIndexSpace<BoundaryType>(
bc_index_type bc_indices = createBoundaryIndexSpace(
exec_space, particles, planes, initial_guess );
return BoundaryCondition<bc_index_type, UserFunctor>(
bc_indices, user_functor, force_update );
}

// Custom index space cases
template <class BCTag, class ExecSpace, class BoundaryParticles>
auto createBoundaryCondition( BCTag, const double value,
BoundaryParticles particles )
{
using memory_space = typename BoundaryParticles::memory_space;
using bc_index_type = BoundaryIndexSpace<memory_space, Custom>;
bc_index_type bc_indices = createBoundaryIndexSpace( particles );
return BoundaryCondition<bc_index_type, BCTag>( value, bc_indices );
}

template <class UserFunctor, class BoundaryParticles>
auto createBoundaryCondition( UserFunctor user_functor,
BoundaryParticles particles,
const bool force_update )
{
using memory_space = typename BoundaryParticles::memory_space;
using bc_index_type = BoundaryIndexSpace<memory_space, Custom>;
bc_index_type bc_indices = createBoundaryIndexSpace( particles );
return BoundaryCondition<bc_index_type, UserFunctor>(
bc_indices, user_functor, force_update );
}

// Wrappers for single plane cases.
template <class BoundaryType, class BCTag, class ExecSpace, class Particles>
auto createBoundaryCondition( BCTag tag, const double value,
streeve marked this conversation as resolved.
Show resolved Hide resolved
ExecSpace exec_space, Particles particles,
BoundaryType plane,
const double initial_guess = 0.5 )
streeve marked this conversation as resolved.
Show resolved Hide resolved
{
std::vector<BoundaryType> plane_vec = { plane };
return createBoundaryCondition( tag, value, exec_space, particles,
plane_vec, initial_guess );
}

template <class UserFunctor, class BoundaryType, class ExecSpace,
class Particles>
auto createBoundaryCondition( UserFunctor user_functor, ExecSpace exec_space,
streeve marked this conversation as resolved.
Show resolved Hide resolved
Particles particles, BoundaryType plane,
const bool force_update,
const double initial_guess = 0.5 )
{
std::vector<BoundaryType> plane_vec = { plane };
return createBoundaryCondition( user_functor, exec_space, particles,
plane_vec, force_update, initial_guess );
}

} // namespace CabanaPD

#endif
Loading