diff --git a/examples/crack_branching.cpp b/examples/crack_branching.cpp index e4c8ddd0..61e346d0 100644 --- a/examples/crack_branching.cpp +++ b/examples/crack_branching.cpp @@ -120,8 +120,8 @@ int main( int argc, char* argv[] ) }; particles->updateParticles( exec_space{}, init_functor ); - // FIXME: use createSolver to switch backend at runtime. - auto cabana_pd = CabanaPD::createSolverFracture( + using iter_tag = Cabana::SerialOpTag; + auto cabana_pd = CabanaPD::createSolverFracture( inputs, particles, force_model, bc, prenotch ); cabana_pd->init_force(); cabana_pd->run(); diff --git a/examples/elastic_wave.cpp b/examples/elastic_wave.cpp index 74271326..852e4918 100644 --- a/examples/elastic_wave.cpp +++ b/examples/elastic_wave.cpp @@ -90,7 +90,8 @@ int main( int argc, char* argv[] ) }; particles->updateParticles( exec_space{}, init_functor ); - auto cabana_pd = CabanaPD::createSolverElastic( + using iter_tag = Cabana::SerialOpTag; + auto cabana_pd = CabanaPD::createSolverElastic( inputs, particles, force_model ); cabana_pd->init_force(); cabana_pd->run(); diff --git a/examples/kalthoff_winkler.cpp b/examples/kalthoff_winkler.cpp index 4c072975..b9bffc09 100644 --- a/examples/kalthoff_winkler.cpp +++ b/examples/kalthoff_winkler.cpp @@ -120,7 +120,9 @@ int main( int argc, char* argv[] ) }; particles->updateParticles( exec_space{}, init_functor ); - auto cabana_pd = CabanaPD::createSolverFracture( + using iter_tag = Cabana::TeamOpTag; + // using iter_tag = Cabana::SerialOpTag; + auto cabana_pd = CabanaPD::createSolverFracture( inputs, particles, force_model, bc, prenotch ); cabana_pd->init_force(); cabana_pd->run(); diff --git a/src/CabanaPD_Force.hpp b/src/CabanaPD_Force.hpp index 13f52459..1a03cd44 100644 --- a/src/CabanaPD_Force.hpp +++ b/src/CabanaPD_Force.hpp @@ -561,10 +561,12 @@ class Force> { } - template + template void computeWeightedVolume( ParticleType& particles, const NeighListType& neigh_list, - const MuView& mu ) + const MuView& mu, + const ParallelType neigh_op_tag ) { auto n_local = particles.n_local; auto x = particles.sliceReferencePosition(); @@ -574,36 +576,33 @@ class Force> Cabana::deep_copy( m, 0.0 ); auto model = _model; - auto weighted_volume = KOKKOS_LAMBDA( const int i ) + auto weighted_volume = KOKKOS_LAMBDA( const int i, const int n ) { - std::size_t num_neighbors = - Cabana::NeighborList::numNeighbor( neigh_list, - i ); - for ( std::size_t n = 0; n < num_neighbors; n++ ) - { - std::size_t j = - Cabana::NeighborList::getNeighbor( - neigh_list, i, n ); + std::size_t j = Cabana::NeighborList::getNeighbor( + neigh_list, i, n ); - // Get the reference positions and displacements. - double xi, r, s; - getDistance( x, u, i, j, xi, r, s ); - // mu is included to account for bond breaking. - double m_j = mu( i, n ) * model.influence_function( xi ) * xi * - xi * vol( j ); - m( i ) += m_j; - } + // Get the reference positions and displacements. + double xi, r, s; + getDistance( x, u, i, j, xi, r, s ); + // mu is included to account for bond breaking. + double m_j = mu( i, n ) * model.influence_function( xi ) * xi * xi * + vol( j ); + m( i ) += m_j; }; Kokkos::RangePolicy policy( 0, n_local ); - Kokkos::parallel_for( "CabanaPD::ForceLPSDamage::computeWeightedVolume", - policy, weighted_volume ); + // Indirect indexing needed to access the broken bonds array. + Cabana::neighbor_parallel_for( + policy, weighted_volume, neigh_list, Cabana::FirstNeighborsTag(), + neigh_op_tag, "CabanaPD::ForceLPSDamage::computeWeightedVolume", + false ); } - template + template void computeDilatation( ParticleType& particles, - const NeighListType& neigh_list, - const MuView& mu ) const + const NeighListType& neigh_list, const MuView& mu, + const ParallelType neigh_op_tag ) const { auto n_local = particles.n_local; const auto x = particles.sliceReferencePosition(); @@ -614,34 +613,30 @@ class Force> auto model = _model; Cabana::deep_copy( theta, 0.0 ); - auto dilatation = KOKKOS_LAMBDA( const int i ) + auto dilatation = KOKKOS_LAMBDA( const int i, const int n ) { - std::size_t num_neighbors = - Cabana::NeighborList::numNeighbor( neigh_list, - i ); - for ( std::size_t n = 0; n < num_neighbors; n++ ) - { - std::size_t j = - Cabana::NeighborList::getNeighbor( - neigh_list, i, n ); + std::size_t j = Cabana::NeighborList::getNeighbor( + neigh_list, i, n ); - // Get the bond distance, displacement, and stretch. - double xi, r, s; - getDistance( x, u, i, j, xi, r, s ); - // mu is included to account for bond breaking. - double theta_i = mu( i, n ) * model.influence_function( xi ) * - s * xi * xi * vol( j ); - // Check if all bonds are broken (m=0) to avoid dividing by - // zero. Alternatively, one could check if this bond mu(i,n) is - // broken, beacuse m=0 only occurs when all bonds are broken. - if ( m( i ) > 0 ) - theta( i ) += 3.0 * theta_i / m( i ); - } + // Get the bond distance, displacement, and stretch. + double xi, r, s; + getDistance( x, u, i, j, xi, r, s ); + // mu is included to account for bond breaking. + double theta_i = mu( i, n ) * model.influence_function( xi ) * s * + xi * xi * vol( j ); + // Check if all bonds are broken (m=0) to avoid dividing by + // zero. Alternatively, one could check if this bond mu(i,n) is + // broken, beacuse m=0 only occurs when all bonds are broken. + if ( m( i ) > 0 ) + theta( i ) += 3.0 * theta_i / m( i ); }; Kokkos::RangePolicy policy( 0, n_local ); - Kokkos::parallel_for( "CabanaPD::ForceLPSDamage::computeDilatation", - policy, dilatation ); + // Indirect indexing needed to access the broken bonds array. + Cabana::neighbor_parallel_for( + policy, dilatation, neigh_list, Cabana::FirstNeighborsTag(), + neigh_op_tag, "CabanaPD::ForceLPSDamage::computeDilatation", + false ); } template > void computeForceFull( ForceType& f, const PosType& x, const PosType& u, const ParticleType& particles, const NeighListType& neigh_list, MuView& mu, - const int n_local, ParallelType& ) const + const int n_local, ParallelType& neigh_op_tag ) const { auto break_coeff = _model.bond_break_coeff; auto theta_coeff = _model.theta_coeff; @@ -660,57 +655,52 @@ class Force> auto theta = particles.sliceDilatation(); auto m = particles.sliceWeightedVolume(); const auto nofail = particles.sliceNoFail(); - auto force_full = KOKKOS_LAMBDA( const int i ) + auto force_full = KOKKOS_LAMBDA( const int i, const int n ) { - std::size_t num_neighbors = - Cabana::NeighborList::numNeighbor( neigh_list, - i ); - for ( std::size_t n = 0; n < num_neighbors; n++ ) - { - double fx_i = 0.0; - double fy_i = 0.0; - double fz_i = 0.0; + double fx_i = 0.0; + double fy_i = 0.0; + double fz_i = 0.0; - std::size_t j = - Cabana::NeighborList::getNeighbor( - neigh_list, i, n ); + std::size_t j = Cabana::NeighborList::getNeighbor( + neigh_list, i, n ); - // Get the reference positions and displacements. - double xi, r, s; - double rx, ry, rz; - getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz ); - - // Break if beyond critical stretch unless in no-fail zone. - if ( r * r >= break_coeff * xi * xi && !nofail( i ) && - !nofail( i ) ) - { - mu( i, n ) = 0; - } - // Check if this bond is broken (mu=0) to ensure m(i) and m(j) - // are both >0 (m=0 only occurs when all bonds are broken) to - // avoid dividing by zero. - else if ( mu( i, n ) > 0 ) - { - const double coeff = - ( theta_coeff * - ( theta( i ) / m( i ) + theta( j ) / m( j ) ) + - s_coeff * s * ( 1.0 / m( i ) + 1.0 / m( j ) ) ) * - model.influence_function( xi ) * xi * vol( j ); - double muij = mu( i, n ); - fx_i = muij * coeff * rx / r; - fy_i = muij * coeff * ry / r; - fz_i = muij * coeff * rz / r; - - f( i, 0 ) += fx_i; - f( i, 1 ) += fy_i; - f( i, 2 ) += fz_i; - } + // Get the reference positions and displacements. + double xi, r, s; + double rx, ry, rz; + getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz ); + + // Break if beyond critical stretch unless in no-fail zone. + if ( r * r >= break_coeff * xi * xi && !nofail( i ) && + !nofail( i ) ) + { + mu( i, n ) = 0; + } + // Check if this bond is broken (mu=0) to ensure m(i) and m(j) + // are both >0 (m=0 only occurs when all bonds are broken) to + // avoid dividing by zero. + else if ( mu( i, n ) > 0 ) + { + const double coeff = + ( theta_coeff * + ( theta( i ) / m( i ) + theta( j ) / m( j ) ) + + s_coeff * s * ( 1.0 / m( i ) + 1.0 / m( j ) ) ) * + model.influence_function( xi ) * xi * vol( j ); + double muij = mu( i, n ); + fx_i = muij * coeff * rx / r; + fy_i = muij * coeff * ry / r; + fz_i = muij * coeff * rz / r; + + f( i, 0 ) += fx_i; + f( i, 1 ) += fy_i; + f( i, 2 ) += fz_i; } }; Kokkos::RangePolicy policy( 0, n_local ); - Kokkos::parallel_for( "CabanaPD::ForceLPSDamage::computeFull", policy, - force_full ); + // Indirect indexing needed to access the broken bonds array. + Cabana::neighbor_parallel_for( + policy, force_full, neigh_list, Cabana::FirstNeighborsTag(), + neigh_op_tag, "CabanaPD::ForceLPSDamage::computeFull", false ); } template > { } + template + void computeWeightedVolume( ParticleType&, const NeighListType&, + const MuView&, const ParallelType ) + { + } + template + void computeDilatation( ParticleType&, const NeighListType&, const MuView&, + const ParallelType ) + { + } + template void computeForceFull( ForceType& f, const PosType& x, const PosType& u, const ParticleType& particles, const NeighListType& neigh_list, MuView& mu, - const int n_local, ParallelType& ) const + const int n_local, ParallelType& neigh_op_tag ) const { auto c = _model.c; auto break_coeff = _model.bond_break_coeff; const auto vol = particles.sliceVolume(); const auto nofail = particles.sliceNoFail(); - auto force_full = KOKKOS_LAMBDA( const int i ) + auto force_full = KOKKOS_LAMBDA( const int i, const int n ) { - std::size_t num_neighbors = - Cabana::NeighborList::numNeighbor( neigh_list, - i ); - for ( std::size_t n = 0; n < num_neighbors; n++ ) - { - double fx_i = 0.0; - double fy_i = 0.0; - double fz_i = 0.0; + double fx_i = 0.0; + double fy_i = 0.0; + double fz_i = 0.0; - std::size_t j = - Cabana::NeighborList::getNeighbor( - neigh_list, i, n ); + std::size_t j = Cabana::NeighborList::getNeighbor( + neigh_list, i, n ); - // Get the reference positions and displacements. - double xi, r, s; - double rx, ry, rz; - getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz ); - - // Break if beyond critical stretch unless in no-fail zone. - if ( r * r >= break_coeff * xi * xi && !nofail( i ) && - !nofail( j ) ) - { - mu( i, n ) = 0; - } - // Else if statement is only for performance. - else if ( mu( i, n ) > 0 ) - { - const double coeff = c * s * vol( j ); - double muij = mu( i, n ); - fx_i = muij * coeff * rx / r; - fy_i = muij * coeff * ry / r; - fz_i = muij * coeff * rz / r; - - f( i, 0 ) += fx_i; - f( i, 1 ) += fy_i; - f( i, 2 ) += fz_i; - } + // Get the reference positions and displacements. + double xi, r, s; + double rx, ry, rz; + getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz ); + + // Break if beyond critical stretch unless in no-fail zone. + if ( r * r >= break_coeff * xi * xi && !nofail( i ) && + !nofail( j ) ) + { + mu( i, n ) = 0; + } + // Else if statement is only for performance. + else if ( mu( i, n ) > 0 ) + { + const double coeff = c * s * vol( j ); + double muij = mu( i, n ); + fx_i = muij * coeff * rx / r; + fy_i = muij * coeff * ry / r; + fz_i = muij * coeff * rz / r; + + f( i, 0 ) += fx_i; + f( i, 1 ) += fy_i; + f( i, 2 ) += fz_i; } }; Kokkos::RangePolicy policy( 0, n_local ); - Kokkos::parallel_for( "CabanaPD::ForcePMBDamage::computeFull", policy, - force_full ); + // Indirect indexing needed to access the broken bonds array. + Cabana::neighbor_parallel_for( + policy, force_full, neigh_list, Cabana::FirstNeighborsTag(), + neigh_op_tag, "CabanaPD::ForcePMBDamage::computeFull", false ); } template +template class SolverElastic { public: @@ -106,7 +106,7 @@ class SolverElastic using neighbor_type = Cabana::VerletList; - using neigh_iter_tag = Cabana::SerialOpTag; + using neigh_iter_tag = IterationTag; using input_type = InputType; SolverElastic( input_type _inputs, @@ -332,14 +332,15 @@ class SolverElastic bool print; }; -template -class SolverFracture - : public SolverElastic +template +class SolverFracture : public SolverElastic { public: - using base_type = - SolverElastic; + using base_type = SolverElastic; using exec_space = typename base_type::exec_space; using memory_space = typename base_type::memory_space; @@ -349,7 +350,7 @@ class SolverFracture using neighbor_type = typename base_type::neighbor_type; using force_model_type = ForceModel; using force_type = typename base_type::force_type; - using neigh_iter_tag = Cabana::SerialOpTag; + using neigh_iter_tag = typename base_type::neigh_iter_tag; using bc_type = BoundaryCondition; using prenotch_type = PrenotchType; using input_type = typename base_type::input_type; @@ -380,10 +381,12 @@ class SolverFracture { init_timer.reset(); // Compute/communicate weighted volume for LPS (does nothing for PMB). - force->computeWeightedVolume( *particles, *neighbors, mu ); + force->computeWeightedVolume( *particles, *neighbors, mu, + neigh_iter_tag{} ); comm->gatherWeightedVolume(); // Compute/communicate dilatation for LPS (does nothing for PMB). - force->computeDilatation( *particles, *neighbors, mu ); + force->computeDilatation( *particles, *neighbors, mu, + neigh_iter_tag{} ); comm->gatherDilatation(); // Compute initial forces. @@ -416,14 +419,16 @@ class SolverFracture // Compute/communicate LPS weighted volume (does nothing for PMB). force_timer.reset(); - force->computeWeightedVolume( *particles, *neighbors, mu ); + force->computeWeightedVolume( *particles, *neighbors, mu, + neigh_iter_tag{} ); force_time += force_timer.seconds(); comm_timer.reset(); comm->gatherWeightedVolume(); comm_time += comm_timer.seconds(); // Compute/communicate LPS dilatation (does nothing for PMB). force_timer.reset(); - force->computeDilatation( *particles, *neighbors, mu ); + force->computeDilatation( *particles, *neighbors, mu, + neigh_iter_tag{} ); force_time += force_timer.seconds(); comm_timer.reset(); comm->gatherDilatation(); @@ -495,27 +500,28 @@ class SolverFracture using base_type::print; }; -template +template auto createSolverElastic( InputsType inputs, std::shared_ptr particles, ForceModel model ) { - return std::make_shared< - SolverElastic>( + return std::make_shared>( inputs, particles, model ); } -template +template auto createSolverFracture( InputsType inputs, std::shared_ptr particles, ForceModel model, BCType bc, PrenotchType prenotch ) { return std::make_shared< - SolverFracture>( inputs, particles, model, bc, - prenotch ); + SolverFracture>( + inputs, particles, model, bc, prenotch ); } /* diff --git a/unit_test/tstForce.hpp b/unit_test/tstForce.hpp index 4aa16454..e621e231 100644 --- a/unit_test/tstForce.hpp +++ b/unit_test/tstForce.hpp @@ -645,12 +645,14 @@ template -void initializeForce( CabanaPD::ForceModel, +template +void initializeForce( CabanaPD::ForceModel, ForceType& force, ParticleType& particles, const NeighborList& neigh_list ) { @@ -659,8 +661,9 @@ void initializeForce( CabanaPD::ForceModel, Kokkos::View mu( "broken_bonds", particles.n_local, max_neighbors ); Kokkos::deep_copy( mu, 1 ); - force.computeWeightedVolume( particles, neigh_list, mu ); - force.computeDilatation( particles, neigh_list, mu ); + force.computeWeightedVolume( particles, neigh_list, mu, + Cabana::SerialOpTag() ); + force.computeDilatation( particles, neigh_list, mu, Cabana::SerialOpTag() ); } template