diff --git a/src/CabanaPD_Comm.hpp b/src/CabanaPD_Comm.hpp index 5ee56eb..d1faad3 100644 --- a/src/CabanaPD_Comm.hpp +++ b/src/CabanaPD_Comm.hpp @@ -195,6 +195,9 @@ struct HaloIds Kokkos::parallel_for( "CabanaPD::Comm::GhostSearch", policy, ghost_search ); Kokkos::fence(); + + // Rebuild if needed. + rebuild( positions ); } template @@ -243,6 +246,7 @@ class Comm int max_export; using memory_space = typename ParticleType::memory_space; + using local_grid_type = typename ParticleType::local_grid_type; using halo_type = Cabana::Halo; using gather_u_type = Cabana::Gather; @@ -262,12 +266,11 @@ class Comm 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 halo_ids( + *local_grid, positions, halo_width, max_export ); // Create the Cabana Halo. halo = std::make_shared( local_grid->globalGrid().comm(), @@ -290,17 +293,6 @@ class Comm auto size() { return mpi_size; } auto rank() { return mpi_rank; } - // Determine which particles should be ghosted, reallocating and recounting - // if needed. - template - auto createHaloIds( const LocalGridType& local_grid, - const PositionSliceType& positions, - const int min_halo_width, const int max_export ) - { - return HaloIds( - 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() @@ -395,6 +387,84 @@ class Comm void gatherTemperature() { gather_temp->apply(); } }; +// Does not inherit because it does not use the same reference halo +// communication pattern. +template +class Comm +{ + public: + using memory_space = typename ParticleType::memory_space; + using local_grid_type = typename ParticleType::local_grid_type; + using halo_type = Cabana::Halo; + std::shared_ptr halo; + HaloIds halo_ids; + + using gather_u_type = + Cabana::Gather; + std::shared_ptr gather_u; + using gather_x_type = + Cabana::Gather; + std::shared_ptr gather_x; + using gather_vol_type = + Cabana::Gather; + std::shared_ptr 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( + *( particles.local_grid ), particles.sliceCurrentPosition(), + halo_width, max_export_guess ) ) + { + auto topology = Cabana::Grid::getTopology( *particles.local_grid ); + halo = std::make_shared( + 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( *halo, particles._aosoa_y ); + gather_x = std::make_shared( + *halo, particles.getReferencePosition().aosoa() ); + gather_vol = + std::make_shared( *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( + 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 diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 84b8732..2fed2d0 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -98,6 +98,8 @@ class Particles 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). @@ -136,9 +138,10 @@ class Particles // FIXME: this is for neighborlist construction. double ghost_mesh_lo[dim]; double ghost_mesh_hi[dim]; - std::shared_ptr< - Cabana::Grid::LocalGrid>> - local_grid; + + using local_grid_type = + Cabana::Grid::LocalGrid>; + std::shared_ptr local_grid; Kokkos::Array dx; int halo_width; @@ -394,7 +397,7 @@ class Particles auto y = Cabana::slice<0>( _aosoa_y, "current_positions" ); auto x = sliceReferencePosition(); auto u = sliceDisplacement(); - Kokkos::RangePolicy policy( 0, n_local + n_ghost ); + Kokkos::RangePolicy policy( 0, size ); auto sum_x_u = KOKKOS_LAMBDA( const std::size_t pid ) { for ( int d = 0; d < 3; d++ ) @@ -405,19 +408,22 @@ class Particles //_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(); }; @@ -466,7 +472,7 @@ class Particles auto time() { return _timer.time(); }; friend class Comm; - friend class Comm; + friend class Comm; protected: aosoa_u_type _aosoa_u; @@ -503,9 +509,11 @@ class Particles 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 @@ -585,12 +593,13 @@ class Particles return Cabana::slice<0>( _aosoa_m, "weighted_volume" ); } - void resize( int new_local, int new_ghost ) + template + void resize( Args&&... args ) { - base_type::resize( new_local, new_ghost ); + base_type::resize( std::forward( 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(); } @@ -638,9 +647,11 @@ class Particles 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 @@ -718,10 +729,11 @@ class Particles return temp_a; } - void resize( int new_local, int new_ghost ) + template + void resize( Args&&... args ) { - base_type::resize( new_local, new_ghost ); - _aosoa_temp.resize( new_local + new_ghost ); + base_type::resize( std::forward( args )... ); + _aosoa_temp.resize( n_reference ); } template @@ -737,6 +749,8 @@ class Particles friend class Comm; friend class Comm; friend class Comm; + friend class Comm; + friend class Comm; protected: void init_temp() @@ -767,6 +781,7 @@ class Particles 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 @@ -832,10 +847,11 @@ class Particles return Cabana::slice<1>( _aosoa_output, "damage" ); } - void resize( int new_local, int new_ghost ) + template + void resize( Args&&... args ) { - base_type::resize( new_local, new_ghost ); - _aosoa_output.resize( new_local + new_ghost ); + base_type::resize( std::forward( args )... ); + _aosoa_output.resize( n_reference ); } template @@ -849,8 +865,10 @@ class Particles friend class Comm; friend class Comm; + friend class Comm; friend class Comm; friend class Comm; + friend class Comm; protected: void init_output() diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 6625086..cceb38d 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -113,6 +113,8 @@ class SolverElastic using heat_transfer_type = HeatTransfer; using contact_type = Force; using contact_model_type = ContactModelType; + using contact_comm_type = + Comm; SolverElastic( input_type _inputs, std::shared_ptr _particles, @@ -150,16 +152,11 @@ class SolverElastic dt = inputs["timestep"]; integrator = std::make_shared( dt ); + std::cout << particles->n_local << " " << particles->n_ghost + << std::endl; // Add ghosts from other MPI ranks. comm = std::make_shared( *particles ); - if constexpr ( is_contact::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 ) @@ -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::value ) + contact_comm = std::make_shared( *particles ); + + std::cout << particles->n_local << " " << particles->n_ghost << " " + << particles->n_contact_ghost << std::endl; + print = print_rank(); if ( print ) { @@ -470,6 +475,7 @@ class SolverElastic // Optional modules. std::shared_ptr heat_transfer; std::shared_ptr contact; + std::shared_ptr contact_comm; contact_model_type contact_model; // Output files. @@ -618,6 +624,11 @@ class SolverFracture // Update ghost particles. comm->gatherDisplacement(); + // Need to ghost current positions for DEM/contact. + if constexpr ( is_contact::value || + is_contact::value ) + contact_comm->gather( *particles ); + // Compute internal forces. updateForce(); @@ -657,6 +668,11 @@ class SolverFracture // Update ghost particles. comm->gatherDisplacement(); + // Need to ghost current positions for DEM/contact. + if constexpr ( is_contact::value || + is_contact::value ) + contact_comm->gather( *particles ); + // Compute internal forces. updateForce(); @@ -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; diff --git a/src/CabanaPD_Types.hpp b/src/CabanaPD_Types.hpp index fa82092..2a1f53e 100644 --- a/src/CabanaPD_Types.hpp +++ b/src/CabanaPD_Types.hpp @@ -24,6 +24,9 @@ struct Fracture // Contact and DEM (contact without PD) tags. +struct Contact +{ +}; struct NoContact { };