Skip to content

Commit

Permalink
Add particle lists for force and position - test with crack branch BC
Browse files Browse the repository at this point in the history
  • Loading branch information
@streeve
streeve committed Jan 8, 2024
1 parent 879961f commit 31f6e81
Show file tree
Hide file tree
Showing 7 changed files with 118 additions and 60 deletions.
17 changes: 15 additions & 2 deletions examples/crack_branching.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,9 +91,22 @@ int main( int argc, char* argv[] )
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 };
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 )
{
// Get a modifiable copy of force.
auto p_f = particles_f.getParticleView( pid );
// Get a copy of the position.
auto p_x = particles_x.getParticle( pid );
auto ypos =
Cabana::get( p_x, CabanaPD::Field::ReferencePosition(), 1 );
auto sign = std::abs( ypos ) / ypos;
Cabana::get( p_f, CabanaPD::Field::Force(), 1 ) += b0 * sign;
};
auto bc =
createBoundaryCondition( CabanaPD::ForceCrackBranchBCTag{},
exec_space{}, *particles, planes, b0 );
createBoundaryCondition( bc_op, exec_space{}, *particles, planes );

auto init_functor = KOKKOS_LAMBDA( const int pid )
{
Expand Down
5 changes: 2 additions & 3 deletions examples/kalthoff_winkler.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,9 +97,8 @@ int main( int argc, char* argv[] )
x_bc - dx, x_bc + dx * 1.25, y_prenotch1 - dx * 0.25,
y_prenotch2 + dx * 0.25, -thickness, thickness );
std::vector<CabanaPD::RegionBoundary> planes = { plane };
auto bc =
createBoundaryCondition( CabanaPD::ForceValueBCTag{}, exec_space{},
*particles, planes, 0.0 );
auto bc = createBoundaryCondition( CabanaPD::ForceValueBCTag{}, 0.0,
exec_space{}, *particles, planes );

auto init_functor = KOKKOS_LAMBDA( const int pid )
{
Expand Down
1 change: 1 addition & 0 deletions src/CabanaPD.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@

#include <CabanaPD_Boundary.hpp>
#include <CabanaPD_Comm.hpp>
#include <CabanaPD_Fields.hpp>
#include <CabanaPD_Force.hpp>
#include <CabanaPD_Input.hpp>
#include <CabanaPD_Integrate.hpp>
Expand Down
87 changes: 45 additions & 42 deletions src/CabanaPD_Boundary.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -132,40 +132,19 @@ struct ForceUpdateBCTag
{
};

struct ForceCrackBranchBCTag
{
};

struct ZeroBCTag
{
};

template <class BCIndexSpace, class BCTag>
struct BoundaryCondition;

template <class BCIndexSpace>
struct BoundaryCondition<BCIndexSpace, ZeroBCTag>
template <class BCIndexSpace, class UserFunctor>
struct BoundaryCondition
{
template <class ExecSpace, class Particles>
void update( ExecSpace, Particles, RegionBoundary )
{
}

template <class ExecSpace, class ParticleType>
void apply( ExecSpace, ParticleType )
{
}
};

template <class BCIndexSpace>
struct BoundaryCondition<BCIndexSpace, ForceValueBCTag>
{
double _value;
BCIndexSpace _index_space;
UserFunctor _user_functor;

BoundaryCondition( const double value, BCIndexSpace bc_index_space )
: _value( value )
, _index_space( bc_index_space )
BoundaryCondition( BCIndexSpace bc_index_space, UserFunctor user )
: _index_space( bc_index_space )
, _user_functor( user )
{
}

Expand All @@ -177,24 +156,34 @@ struct BoundaryCondition<BCIndexSpace, ForceValueBCTag>
}

template <class ExecSpace, class ParticleType>
void apply( ExecSpace, ParticleType& particles )
void apply( ExecSpace, ParticleType& )
{
auto f = particles.sliceForce();
auto x = particles.sliceReferencePosition();
auto index_space = _index_space._view;
Kokkos::RangePolicy<ExecSpace> policy( 0, index_space.size() );
auto value = _value;
Kokkos::parallel_for(
"CabanaPD::BC::apply", policy, KOKKOS_LAMBDA( const int b ) {
auto pid = index_space( b );
for ( int d = 0; d < 3; d++ )
f( pid, d ) = value;
_user_functor( pid );
} );
}
};

template <class BCIndexSpace>
struct BoundaryCondition<BCIndexSpace, ForceUpdateBCTag>
struct BoundaryCondition<BCIndexSpace, ZeroBCTag>
{
template <class ExecSpace, class Particles>
void update( ExecSpace, Particles, RegionBoundary )
{
}

template <class ExecSpace, class ParticleType>
void apply( ExecSpace, ParticleType )
{
}
};

template <class BCIndexSpace>
struct BoundaryCondition<BCIndexSpace, ForceValueBCTag>
{
double _value;
BCIndexSpace _index_space;
Expand Down Expand Up @@ -223,13 +212,13 @@ struct BoundaryCondition<BCIndexSpace, ForceUpdateBCTag>
"CabanaPD::BC::apply", policy, KOKKOS_LAMBDA( const int b ) {
auto pid = index_space( b );
for ( int d = 0; d < 3; d++ )
f( pid, d ) += value;
f( pid, d ) = value;
} );
}
};

template <class BCIndexSpace>
struct BoundaryCondition<BCIndexSpace, ForceCrackBranchBCTag>
struct BoundaryCondition<BCIndexSpace, ForceUpdateBCTag>
{
double _value;
BCIndexSpace _index_space;
Expand All @@ -251,25 +240,23 @@ struct BoundaryCondition<BCIndexSpace, ForceCrackBranchBCTag>
void apply( ExecSpace, ParticleType& particles )
{
auto f = particles.sliceForce();
auto x = particles.sliceReferencePosition();
auto index_space = _index_space._view;
Kokkos::RangePolicy<ExecSpace> policy( 0, index_space.size() );
auto value = _value;
Kokkos::parallel_for(
"CabanaPD::BC::apply", policy, KOKKOS_LAMBDA( const int b ) {
auto pid = index_space( b );
// This is specifically for the crack branching.
auto sign = std::abs( x( pid, 1 ) ) / x( pid, 1 );
f( pid, 1 ) += value * sign;
for ( int d = 0; d < 3; d++ )
f( pid, d ) += value;
} );
}
};

// FIXME: relatively large initial guess for allocation.
template <class BoundaryType, class BCTag, class ExecSpace, class Particles>
auto createBoundaryCondition( BCTag, ExecSpace exec_space, Particles particles,
auto createBoundaryCondition( BCTag, const double value, ExecSpace exec_space,
Particles particles,
std::vector<BoundaryType> planes,
const double value,
const double initial_guess = 0.5 )
{
using memory_space = typename Particles::memory_space;
Expand All @@ -279,6 +266,22 @@ auto createBoundaryCondition( BCTag, ExecSpace exec_space, Particles particles,
return BoundaryCondition<bc_index_type, BCTag>( value, bc_indices );
}

// FIXME: relatively large initial guess for allocation.
template <class UserFunctor, class BoundaryType, class ExecSpace,
class Particles>
auto createBoundaryCondition( UserFunctor user_functor, ExecSpace exec_space,
Particles particles,
std::vector<BoundaryType> planes,
const double initial_guess = 0.5 )
{
using memory_space = typename Particles::memory_space;
using bc_index_type = BoundaryIndexSpace<memory_space, BoundaryType>;
bc_index_type bc_indices = createBoundaryIndexSpace<BoundaryType>(
exec_space, particles, planes, initial_guess );
return BoundaryCondition<bc_index_type, UserFunctor>( bc_indices,
user_functor );
}

template <class MemorySpace>
auto createBoundaryCondition( ZeroBCTag )
{
Expand Down
3 changes: 1 addition & 2 deletions src/CabanaPD_Comm.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -276,9 +276,8 @@ class Comm<ParticleType, PMB>

// Only use this interface because we don't need to recommunicate
// positions, volumes, or no-fail region.
Cabana::gather( *halo, particles._aosoa_x );
Cabana::gather( *halo, particles.getReferencePosition().aosoa() );
Cabana::gather( *halo, particles._aosoa_vol );
Cabana::gather( *halo, particles._aosoa_nofail );

gather_u = std::make_shared<gather_u_type>( *halo, particles._aosoa_u );
gather_u->apply();
Expand Down
32 changes: 32 additions & 0 deletions src/CabanaPD_Fields.hpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,32 @@
#ifndef FIELDS_HPP
#define FIELDS_HPP

#include <Cabana_Core.hpp>

namespace CabanaPD
{
//---------------------------------------------------------------------------//
// Fields.
//---------------------------------------------------------------------------//
namespace Field
{

struct ReferencePosition : Cabana::Field::Position<3>
{
static std::string label() { return "reference_positions"; }
};

struct Force : Cabana::Field::Vector<double, 3>
{
static std::string label() { return "forces"; }
};

struct NoFail : Cabana::Field::Scalar<int>
{
static std::string label() { return "no_fail"; }
};

} // namespace Field
} // namespace CabanaPD

#endif
33 changes: 22 additions & 11 deletions src/CabanaPD_Particles.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,13 +106,16 @@ class Particles<MemorySpace, PMB, Dimension>

// FIXME: add vector length.
// FIXME: enable variable aosoa.
using aosoa_x_type = Cabana::AoSoA<vector_type, memory_space, 1>;
using aosoa_u_type = Cabana::AoSoA<vector_type, memory_space, 1>;
using aosoa_y_type = Cabana::AoSoA<vector_type, memory_space, 1>;
using aosoa_f_type = Cabana::AoSoA<vector_type, memory_space, 1>;
using aosoa_vol_type = Cabana::AoSoA<scalar_type, memory_space, 1>;
using aosoa_nofail_type = Cabana::AoSoA<int_type, memory_space, 1>;
using aosoa_other_type = Cabana::AoSoA<other_types, memory_space>;
using plist_x_type =
Cabana::ParticleList<memory_space, 1,
CabanaPD::Field::ReferencePosition>;
using plist_f_type =
Cabana::ParticleList<memory_space, 1, CabanaPD::Field::Force>;

// Per type.
int n_types = 1;
Expand Down Expand Up @@ -153,6 +156,8 @@ class Particles<MemorySpace, PMB, Dimension>
std::array<double, dim> high_corner,
const std::array<int, dim> num_cells, const int max_halo_width )
: halo_width( max_halo_width )
, _plist_x( "positions" )
, _plist_f( "forces" )
{
createDomain( low_corner, high_corner, num_cells );
createParticles( exec_space );
Expand Down Expand Up @@ -205,6 +210,8 @@ class Particles<MemorySpace, PMB, Dimension>

int particles_per_cell = 1;
int num_particles = particles_per_cell * owned_cells.size();

// Use default aosoa construction and resize.
resize( num_particles, 0 );

auto x = sliceReferencePosition();
Expand Down Expand Up @@ -275,7 +282,7 @@ class Particles<MemorySpace, PMB, Dimension>
local_num_create );
n_local = local_num_create;
resize( n_local, 0 );
size = _aosoa_x.size();
size = _plist_x.size();

// Not using Allreduce because global count is only used for printing.
MPI_Reduce( &n_local, &n_global, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, 0,
Expand All @@ -293,11 +300,11 @@ class Particles<MemorySpace, PMB, Dimension>

auto sliceReferencePosition()
{
return Cabana::slice<0>( _aosoa_x, "reference_positions" );
return _plist_x.slice( CabanaPD::Field::ReferencePosition() );
}
auto sliceReferencePosition() const
{
return Cabana::slice<0>( _aosoa_x, "reference_positions" );
return _plist_x.slice( CabanaPD::Field::ReferencePosition() );
}
auto sliceCurrentPosition()
{
Expand All @@ -319,7 +326,7 @@ class Particles<MemorySpace, PMB, Dimension>
{
return Cabana::slice<0>( _aosoa_u, "displacements" );
}
auto sliceForce() { return Cabana::slice<0>( _aosoa_f, "forces" ); }
auto sliceForce() { return _plist_f.slice( CabanaPD::Field::Force() ); }
auto sliceForceAtomic()
{
auto f = sliceForce();
Expand Down Expand Up @@ -370,6 +377,9 @@ class Particles<MemorySpace, PMB, Dimension>
return Cabana::slice<0>( _aosoa_nofail, "no_fail_region" );
}

auto getForce() { return _plist_f; }
auto getReferencePosition() { return _plist_x; }

void updateCurrentPosition()
{
// Not using slice function because this is called inside.
Expand All @@ -391,14 +401,14 @@ class Particles<MemorySpace, PMB, Dimension>
n_local = new_local;
n_ghost = new_ghost;

_aosoa_x.resize( new_local + new_ghost );
_plist_x.aosoa().resize( new_local + new_ghost );
_aosoa_u.resize( new_local + new_ghost );
_aosoa_y.resize( new_local + new_ghost );
_aosoa_vol.resize( new_local + new_ghost );
_aosoa_f.resize( new_local );
_plist_f.aosoa().resize( new_local );
_aosoa_other.resize( new_local );
_aosoa_nofail.resize( new_local + new_ghost );
size = _aosoa_x.size();
size = _plist_x.size();
};

auto getPosition( const bool use_reference )
Expand Down Expand Up @@ -433,14 +443,15 @@ class Particles<MemorySpace, PMB, Dimension>
friend class Comm<self_type, PMB>;

protected:
aosoa_x_type _aosoa_x;
aosoa_u_type _aosoa_u;
aosoa_y_type _aosoa_y;
aosoa_f_type _aosoa_f;
aosoa_vol_type _aosoa_vol;
aosoa_nofail_type _aosoa_nofail;
aosoa_other_type _aosoa_other;

plist_x_type _plist_x;
plist_f_type _plist_f;

#ifdef Cabana_ENABLE_HDF5
Cabana::Experimental::HDF5ParticleOutput::HDF5Config h5_config;
#endif
Expand Down

0 comments on commit 31f6e81

Please sign in to comment.