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

Add contact ghosts #157

Draft
wants to merge 3 commits into
base: main
Choose a base branch
from
Draft

Add contact ghosts #157

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
8c70864
Add contact ghosts
streeve Dec 6, 2024
542014f
fixup: contact comm compile issues
streeve Dec 6, 2024
f5f8ddf
fixup: communicate x, u, and volume for contact
streeve Dec 6, 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
102 changes: 86 additions & 16 deletions src/CabanaPD_Comm.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -195,6 +195,9 @@ struct HaloIds
Kokkos::parallel_for( "CabanaPD::Comm::GhostSearch", policy,
ghost_search );
Kokkos::fence();

// Rebuild if needed.
rebuild( positions );
}

template <class PositionSliceType>
Expand Down Expand Up @@ -243,6 +246,7 @@ class Comm<ParticleType, PMB, TemperatureIndependent>
int max_export;

using memory_space = typename ParticleType::memory_space;
using local_grid_type = typename ParticleType::local_grid_type;
using halo_type = Cabana::Halo<memory_space>;
using gather_u_type =
Cabana::Gather<halo_type, typename ParticleType::aosoa_u_type>;
Expand All @@ -262,12 +266,11 @@ class Comm<ParticleType, PMB, TemperatureIndependent>
auto halo_width = local_grid->haloCellWidth();
auto topology = Cabana::Grid::getTopology( *local_grid );

// Determine which particles need to be ghosted to neighbors.
// Determine which particles should be ghosted, reallocating and
// recounting if needed.
// FIXME: set halo width based on cutoff distance.
auto halo_ids =
createHaloIds( *local_grid, positions, halo_width, max_export );
// Rebuild if needed.
halo_ids.rebuild( positions );
HaloIds<memory_space, local_grid_type> halo_ids(
*local_grid, positions, halo_width, max_export );

// Create the Cabana Halo.
halo = std::make_shared<halo_type>( local_grid->globalGrid().comm(),
Expand All @@ -290,17 +293,6 @@ class Comm<ParticleType, PMB, TemperatureIndependent>
auto size() { return mpi_size; }
auto rank() { return mpi_rank; }

// Determine which particles should be ghosted, reallocating and recounting
// if needed.
template <class LocalGridType, class PositionSliceType>
auto createHaloIds( const LocalGridType& local_grid,
const PositionSliceType& positions,
const int min_halo_width, const int max_export )
{
return HaloIds<typename PositionSliceType::memory_space, LocalGridType>(
local_grid, positions, min_halo_width, max_export );
}

// We assume here that the particle count has not changed and no resize
// is necessary.
void gatherDisplacement()
Expand Down Expand Up @@ -395,6 +387,84 @@ class Comm<ParticleType, PMB, TemperatureDependent>
void gatherTemperature() { gather_temp->apply(); }
};

// Does not inherit because it does not use the same reference halo
// communication pattern.
template <class ParticleType>
class Comm<ParticleType, Contact, TemperatureIndependent>
{
public:
using memory_space = typename ParticleType::memory_space;
using local_grid_type = typename ParticleType::local_grid_type;
using halo_type = Cabana::Halo<memory_space>;
std::shared_ptr<halo_type> halo;
HaloIds<memory_space, local_grid_type> halo_ids;

using gather_u_type =
Cabana::Gather<halo_type, typename ParticleType::aosoa_u_type>;
std::shared_ptr<gather_u_type> gather_u;
using gather_x_type =
Cabana::Gather<halo_type,
typename ParticleType::plist_x_type::aosoa_type>;
std::shared_ptr<gather_x_type> gather_x;
using gather_vol_type =
Cabana::Gather<halo_type, typename ParticleType::aosoa_vol_type>;
std::shared_ptr<gather_vol_type> gather_vol;

// Note this initial guess is small because this is often used for very
// short range interactions.
Comm( ParticleType& particles, int halo_width = 1,
int max_export_guess = 10 )
: halo_ids( HaloIds<memory_space, local_grid_type>(
*( particles.local_grid ), particles.sliceCurrentPosition(),
halo_width, max_export_guess ) )
{
auto topology = Cabana::Grid::getTopology( *particles.local_grid );
halo = std::make_shared<halo_type>(
particles.local_grid->globalGrid().comm(), particles.n_reference,
halo_ids._ids, halo_ids._destinations, topology );
std::cout << "halo " << halo->numLocal() << " " << halo->numGhost()
<< std::endl;
// We use n_ghost here as the "local" halo count because these current
// frame ghosts are built on top of the existing, static, reference
// frame ghosts.
particles.resize( particles.n_local, particles.n_ghost,
halo->numGhost() );

gather_u = std::make_shared<gather_u_type>( *halo, particles._aosoa_y );
gather_x = std::make_shared<gather_x_type>(
*halo, particles.getReferencePosition().aosoa() );
gather_vol =
std::make_shared<gather_vol_type>( *halo, particles._aosoa_vol );
}

// This is a dynamic gather step where the steering vector needs to be
// recomputed.
void gather( ParticleType& particles )
{
// Get the current position. Note this is necessary to get the up to
// date current position.
auto y = particles.sliceCurrentPosition();
// Determine which particles need to be ghosted to neighbors for the
// current positions.
halo_ids.build( y );

auto topology = Cabana::Grid::getTopology( *particles.local_grid );
// FIXME: missing a build() interface
halo = std::make_shared<halo_type>(
particles.local_grid->globalGrid().comm(), particles.n_reference,
halo_ids._ids, halo_ids._destinations, topology );
particles.resize( particles.n_local, particles.n_ghost,
halo->numGhost() );

gather_u->reserve( *halo, particles._aosoa_u );
gather_u->apply();
gather_x->reserve( *halo, particles.getReferencePosition().aosoa() );
gather_x->apply();
gather_vol->reserve( *halo, particles._aosoa_vol );
gather_vol->apply();
}
};

} // namespace CabanaPD

#endif
66 changes: 42 additions & 24 deletions src/CabanaPD_Particles.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -98,6 +98,8 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, BaseOutput, Dimension>
unsigned long long int n_global = 0;
std::size_t n_local = 0;
std::size_t n_ghost = 0;
std::size_t n_reference = 0;
std::size_t n_contact_ghost = 0;
std::size_t size = 0;

// x, u, f (vector matching system dimension).
Expand Down Expand Up @@ -136,9 +138,10 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, BaseOutput, Dimension>
// FIXME: this is for neighborlist construction.
double ghost_mesh_lo[dim];
double ghost_mesh_hi[dim];
std::shared_ptr<
Cabana::Grid::LocalGrid<Cabana::Grid::UniformMesh<double, dim>>>
local_grid;

using local_grid_type =
Cabana::Grid::LocalGrid<Cabana::Grid::UniformMesh<double, dim>>;
std::shared_ptr<local_grid_type> local_grid;
Kokkos::Array<double, dim> dx;

int halo_width;
Expand Down Expand Up @@ -394,7 +397,7 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, BaseOutput, Dimension>
auto y = Cabana::slice<0>( _aosoa_y, "current_positions" );
auto x = sliceReferencePosition();
auto u = sliceDisplacement();
Kokkos::RangePolicy<execution_space> policy( 0, n_local + n_ghost );
Kokkos::RangePolicy<execution_space> policy( 0, size );
auto sum_x_u = KOKKOS_LAMBDA( const std::size_t pid )
{
for ( int d = 0; d < 3; d++ )
Expand All @@ -405,19 +408,22 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, BaseOutput, Dimension>
//_timer.stop();
}

void resize( int new_local, int new_ghost )
void resize( int new_local, int new_ghost, int new_contact_ghost = 0 )
{
_timer.start();
n_local = new_local;
n_ghost = 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 );
_plist_f.aosoa().resize( new_local );
_aosoa_other.resize( new_local );
_aosoa_nofail.resize( new_local + new_ghost );
n_reference = n_local + n_ghost;
n_contact_ghost = new_contact_ghost;
int n_all = n_reference + n_contact_ghost;

_plist_x.aosoa().resize( n_all );
_aosoa_u.resize( n_all );
_aosoa_y.resize( n_all );
_aosoa_vol.resize( n_all );
_plist_f.aosoa().resize( n_local );
_aosoa_other.resize( n_local );
_aosoa_nofail.resize( n_reference );
size = _plist_x.size();
_timer.stop();
};
Expand Down Expand Up @@ -466,7 +472,7 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, BaseOutput, Dimension>
auto time() { return _timer.time(); };

friend class Comm<self_type, PMB, TemperatureIndependent>;
friend class Comm<self_type, PMB, TemperatureDependent>;
friend class Comm<self_type, Contact, TemperatureIndependent>;

protected:
aosoa_u_type _aosoa_u;
Expand Down Expand Up @@ -503,9 +509,11 @@ class Particles<MemorySpace, LPS, TemperatureIndependent, BaseOutput, Dimension>
using base_type::dim;

// Per particle.
using base_type::n_contact_ghost;
using base_type::n_ghost;
using base_type::n_global;
using base_type::n_local;
using base_type::n_reference;
using base_type::size;

// These are split since weighted volume only needs to be communicated once
Expand Down Expand Up @@ -585,12 +593,13 @@ class Particles<MemorySpace, LPS, TemperatureIndependent, BaseOutput, Dimension>
return Cabana::slice<0>( _aosoa_m, "weighted_volume" );
}

void resize( int new_local, int new_ghost )
template <typename... Args>
void resize( Args&&... args )
{
base_type::resize( new_local, new_ghost );
base_type::resize( std::forward<Args>( args )... );
_timer.start();
_aosoa_theta.resize( new_local + new_ghost );
_aosoa_m.resize( new_local + new_ghost );
_aosoa_theta.resize( n_reference );
_aosoa_m.resize( n_reference );
_timer.stop();
}

Expand Down Expand Up @@ -638,9 +647,11 @@ class Particles<MemorySpace, PMB, TemperatureDependent, BaseOutput, Dimension>
using base_type::dim;

// Per particle.
using base_type::n_contact_ghost;
using base_type::n_ghost;
using base_type::n_global;
using base_type::n_local;
using base_type::n_reference;
using base_type::size;

// These are split since weighted volume only needs to be communicated once
Expand Down Expand Up @@ -718,10 +729,11 @@ class Particles<MemorySpace, PMB, TemperatureDependent, BaseOutput, Dimension>
return temp_a;
}

void resize( int new_local, int new_ghost )
template <typename... Args>
void resize( Args&&... args )
{
base_type::resize( new_local, new_ghost );
_aosoa_temp.resize( new_local + new_ghost );
base_type::resize( std::forward<Args>( args )... );
_aosoa_temp.resize( n_reference );
}

template <typename... OtherFields>
Expand All @@ -737,6 +749,8 @@ class Particles<MemorySpace, PMB, TemperatureDependent, BaseOutput, Dimension>
friend class Comm<self_type, LPS, TemperatureIndependent>;
friend class Comm<self_type, PMB, TemperatureDependent>;
friend class Comm<self_type, LPS, TemperatureDependent>;
friend class Comm<self_type, Contact, TemperatureIndependent>;
friend class Comm<self_type, Contact, TemperatureDependent>;

protected:
void init_temp()
Expand Down Expand Up @@ -767,6 +781,7 @@ class Particles<MemorySpace, ModelType, ThermalType, EnergyOutput, Dimension>
using base_type::n_ghost;
using base_type::n_global;
using base_type::n_local;
using base_type::n_reference;
using base_type::size;

// energy, damage
Expand Down Expand Up @@ -832,10 +847,11 @@ class Particles<MemorySpace, ModelType, ThermalType, EnergyOutput, Dimension>
return Cabana::slice<1>( _aosoa_output, "damage" );
}

void resize( int new_local, int new_ghost )
template <typename... Args>
void resize( Args&&... args )
{
base_type::resize( new_local, new_ghost );
_aosoa_output.resize( new_local + new_ghost );
base_type::resize( std::forward<Args>( args )... );
_aosoa_output.resize( n_reference );
}

template <typename... OtherFields>
Expand All @@ -849,8 +865,10 @@ class Particles<MemorySpace, ModelType, ThermalType, EnergyOutput, Dimension>

friend class Comm<self_type, PMB, TemperatureIndependent>;
friend class Comm<self_type, LPS, TemperatureIndependent>;
friend class Comm<self_type, Contact, TemperatureIndependent>;
friend class Comm<self_type, PMB, TemperatureDependent>;
friend class Comm<self_type, LPS, TemperatureDependent>;
friend class Comm<self_type, Contact, TemperatureDependent>;

protected:
void init_output()
Expand Down
31 changes: 24 additions & 7 deletions src/CabanaPD_Solver.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -113,6 +113,8 @@ class SolverElastic
using heat_transfer_type = HeatTransfer<memory_space, force_model_type>;
using contact_type = Force<memory_space, ContactModelType>;
using contact_model_type = ContactModelType;
using contact_comm_type =
Comm<particle_type, Contact, TemperatureIndependent>;

SolverElastic( input_type _inputs,
std::shared_ptr<particle_type> _particles,
Expand Down Expand Up @@ -150,16 +152,11 @@ class SolverElastic
dt = inputs["timestep"];
integrator = std::make_shared<integrator_type>( dt );

std::cout << particles->n_local << " " << particles->n_ghost
<< std::endl;
// Add ghosts from other MPI ranks.
comm = std::make_shared<comm_type>( *particles );

if constexpr ( is_contact<contact_model_type>::value )
{
if ( comm->size() > 1 )
throw std::runtime_error(
"Contact with MPI is currently disabled." );
}

// Update temperature ghost size if needed.
if constexpr ( is_temperature_dependent<
typename force_model_type::thermal_type>::value )
Expand Down Expand Up @@ -194,6 +191,14 @@ class SolverElastic
inputs["half_neigh"], force->_neigh_list, force_model );
}

// Purposely delay creating initial contact ghosts until after reference
// neighbor list.
if constexpr ( is_contact<contact_model_type>::value )
contact_comm = std::make_shared<contact_comm_type>( *particles );

std::cout << particles->n_local << " " << particles->n_ghost << " "
<< particles->n_contact_ghost << std::endl;

print = print_rank();
if ( print )
{
Expand Down Expand Up @@ -470,6 +475,7 @@ class SolverElastic
// Optional modules.
std::shared_ptr<heat_transfer_type> heat_transfer;
std::shared_ptr<contact_type> contact;
std::shared_ptr<contact_comm_type> contact_comm;
contact_model_type contact_model;

// Output files.
Expand Down Expand Up @@ -618,6 +624,11 @@ class SolverFracture
// Update ghost particles.
comm->gatherDisplacement();

// Need to ghost current positions for DEM/contact.
if constexpr ( is_contact<contact_model_type>::value ||
is_contact<force_model_type>::value )
contact_comm->gather( *particles );

// Compute internal forces.
updateForce();

Expand Down Expand Up @@ -657,6 +668,11 @@ class SolverFracture
// Update ghost particles.
comm->gatherDisplacement();

// Need to ghost current positions for DEM/contact.
if constexpr ( is_contact<contact_model_type>::value ||
is_contact<force_model_type>::value )
contact_comm->gather( *particles );

// Compute internal forces.
updateForce();

Expand Down Expand Up @@ -716,6 +732,7 @@ class SolverFracture
protected:
using base_type::comm;
using base_type::contact;
using base_type::contact_comm;
using base_type::force;
using base_type::inputs;
using base_type::integrator;
Expand Down
Loading
Loading