diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 080ec4b4..c922bba0 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,8 +1,10 @@ configure_file(CabanaPD_config.hpp.cmakein CabanaPD_config.hpp @ONLY) file(GLOB HEADERS GLOB *.hpp) +file(GLOB FORCE_HEADERS GLOB force/*.hpp) install(FILES ${HEADERS} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) +install(FILES ${FORCE_HEADERS} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/force) install(FILES ${CMAKE_CURRENT_BINARY_DIR}/CabanaPD_config.hpp DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) add_library(CabanaPD INTERFACE) diff --git a/src/CabanaPD.hpp b/src/CabanaPD.hpp index bd094231..b792d878 100644 --- a/src/CabanaPD.hpp +++ b/src/CabanaPD.hpp @@ -17,6 +17,7 @@ #include #include #include +#include #include #include #include @@ -26,4 +27,9 @@ #include #include +#include +#include +#include +#include + #endif diff --git a/src/CabanaPD_Force.hpp b/src/CabanaPD_Force.hpp index fce2c876..3390f58e 100644 --- a/src/CabanaPD_Force.hpp +++ b/src/CabanaPD_Force.hpp @@ -62,260 +62,10 @@ #include -#include -#include #include -#include namespace CabanaPD { -/****************************************************************************** - Force models -******************************************************************************/ -struct BaseForceModel -{ - double delta; - - BaseForceModel(){}; - BaseForceModel( const double _delta ) - : delta( _delta ){}; - - BaseForceModel( const BaseForceModel& model ) { delta = model.delta; } -}; - -template -struct ForceModel; - -/* LPS */ -template <> -struct ForceModel : public BaseForceModel -{ - using base_type = BaseForceModel; - using base_model = LPS; - using fracture_type = Elastic; - - using BaseForceModel::delta; - - int influence_type; - - double K; - double G; - double theta_coeff; - double s_coeff; - - ForceModel(){}; - ForceModel( const double _delta, const double _K, const double _G, - const int _influence = 0 ) - : base_type( _delta ) - , influence_type( _influence ) - { - set_param( _delta, _K, _G ); - } - - void set_param( const double _delta, const double _K, const double _G ) - { - delta = _delta; - K = _K; - G = _G; - - theta_coeff = 3.0 * K - 5.0 * G; - s_coeff = 15.0 * G; - } - - KOKKOS_INLINE_FUNCTION double influence_function( double xi ) const - { - if ( influence_type == 1 ) - return 1.0 / xi; - else - return 1.0; - } -}; - -template <> -struct ForceModel : public ForceModel -{ - using base_type = ForceModel; - using base_model = LPS; - using fracture_type = Fracture; - - using base_type::delta; - using base_type::G; - using base_type::influence_type; - using base_type::K; - using base_type::s_coeff; - using base_type::theta_coeff; - double G0; - double s0; - double bond_break_coeff; - - ForceModel() {} - ForceModel( const double _delta, const double _K, const double _G, - const double _G0, const int _influence = 0 ) - : base_type( _delta, _K, _G, _influence ) - { - set_param( _delta, _K, _G, _G0 ); - } - - void set_param( const double _delta, const double _K, const double _G, - const double _G0 ) - { - base_type::set_param( _delta, _K, _G ); - G0 = _G0; - if ( influence_type == 1 ) - { - s0 = sqrt( 5.0 * G0 / 9.0 / K / delta ); // 1/xi - } - else - { - s0 = sqrt( 8.0 * G0 / 15.0 / K / delta ); // 1 - } - bond_break_coeff = ( 1.0 + s0 ) * ( 1.0 + s0 ); - } -}; - -template <> -struct ForceModel : public ForceModel -{ - using base_type = ForceModel; - using base_model = LPS; - using fracture_type = Elastic; - - using base_type::base_type; - - using base_type::delta; - using base_type::G; - using base_type::influence_type; - using base_type::K; - using base_type::s_coeff; - using base_type::theta_coeff; -}; - -template <> -struct ForceModel : public ForceModel -{ - using base_type = ForceModel; - using base_model = LPS; - using fracture_type = Fracture; - - using base_type::base_type; - - using base_type::delta; - using base_type::G; - using base_type::influence_type; - using base_type::K; - using base_type::s_coeff; - using base_type::theta_coeff; - - using base_type::bond_break_coeff; - using base_type::G0; - using base_type::s0; -}; - -/* PMB */ -template <> -struct ForceModel : public BaseForceModel -{ - using base_type = BaseForceModel; - using base_model = PMB; - using fracture_type = Elastic; - - using BaseForceModel::delta; - - double c; - double K; - - ForceModel(){}; - ForceModel( const double delta, const double K ) - : base_type( delta ) - { - set_param( delta, K ); - } - - ForceModel( const ForceModel& model ) - : base_type( model ) - { - c = model.c; - K = model.K; - } - - void set_param( const double _delta, const double _K ) - { - delta = _delta; - K = _K; - c = 18.0 * K / ( 3.141592653589793 * delta * delta * delta * delta ); - } -}; - -template <> -struct ForceModel : public ForceModel -{ - using base_type = ForceModel; - using base_model = PMB; - using fracture_type = Fracture; - - using base_type::c; - using base_type::delta; - using base_type::K; - double G0; - double s0; - double bond_break_coeff; - - ForceModel() {} - ForceModel( const double delta, const double K, const double G0 ) - : base_type( delta, K ) - { - set_param( delta, K, G0 ); - } - - ForceModel( const ForceModel& model ) - : base_type( model ) - { - G0 = model.G0; - s0 = model.s0; - bond_break_coeff = model.bond_break_coeff; - } - - void set_param( const double _delta, const double _K, const double _G0 ) - { - base_type::set_param( _delta, _K ); - G0 = _G0; - s0 = sqrt( 5.0 * G0 / 9.0 / K / delta ); - bond_break_coeff = ( 1.0 + s0 ) * ( 1.0 + s0 ); - } -}; - -template <> -struct ForceModel : public ForceModel -{ - using base_type = ForceModel; - using base_model = PMB; - using fracture_type = Elastic; - - using base_type::base_type; - - using base_type::c; - using base_type::delta; - using base_type::K; -}; - -template <> -struct ForceModel : public ForceModel -{ - using base_type = ForceModel; - using base_model = PMB; - using fracture_type = Fracture; - - using base_type::base_type; - - using base_type::c; - using base_type::delta; - using base_type::K; - - using base_type::bond_break_coeff; - using base_type::G0; - using base_type::s0; -}; - /****************************************************************************** Force helper functions. ******************************************************************************/ @@ -374,835 +124,13 @@ getLinearizedDistance( const PosType& x, const PosType& u, const int i, getLinearizedDistanceComponents( x, u, i, j, xi, s, xi_x, xi_y, xi_z ); } -/****************************************************************************** - Force computation LPS -******************************************************************************/ +// Forward declaration. template class Force; -template -class Force> -{ - protected: - bool _half_neigh; - ForceModel _model; - - public: - using exec_space = ExecutionSpace; - - Force( const bool half_neigh, const ForceModel model ) - : _half_neigh( half_neigh ) - , _model( model ) - { - } - - template - void computeWeightedVolume( ParticleType& particles, - const NeighListType& neigh_list, - const ParallelType neigh_op_tag ) - { - auto n_local = particles.n_local; - auto x = particles.sliceReferencePosition(); - auto u = particles.sliceDisplacement(); - const auto vol = particles.sliceVolume(); - auto m = particles.sliceWeightedVolume(); - Cabana::deep_copy( m, 0.0 ); - auto model = _model; - - auto weighted_volume = KOKKOS_LAMBDA( const int i, const int j ) - { - // Get the reference positions and displacements. - double xi, r, s; - getDistance( x, u, i, j, xi, r, s ); - double m_j = model.influence_function( xi ) * xi * xi * vol( j ); - m( i ) += m_j; - }; - - Kokkos::RangePolicy policy( 0, n_local ); - Cabana::neighbor_parallel_for( - policy, weighted_volume, neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::ForceLPS::computeWeightedVolume" ); - } - - template - void computeDilatation( ParticleType& particles, - const NeighListType& neigh_list, - const ParallelType neigh_op_tag ) const - { - auto n_local = particles.n_local; - const auto x = particles.sliceReferencePosition(); - auto u = particles.sliceDisplacement(); - const auto vol = particles.sliceVolume(); - auto m = particles.sliceWeightedVolume(); - auto theta = particles.sliceDilatation(); - auto model = _model; - Cabana::deep_copy( theta, 0.0 ); - - auto dilatation = KOKKOS_LAMBDA( const int i, const int j ) - { - // Get the bond distance, displacement, and stretch. - double xi, r, s; - getDistance( x, u, i, j, xi, r, s ); - double theta_i = - model.influence_function( xi ) * s * xi * xi * vol( j ); - theta( i ) += 3.0 * theta_i / m( i ); - }; - - Kokkos::RangePolicy policy( 0, n_local ); - Cabana::neighbor_parallel_for( - policy, dilatation, neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::ForceLPS::computeDilatation" ); - } - - template - void computeForceFull( ForceType& f, const PosType& x, const PosType& u, - const ParticleType& particles, - const NeighListType& neigh_list, const int n_local, - ParallelType& neigh_op_tag ) const - { - auto theta_coeff = _model.theta_coeff; - auto s_coeff = _model.s_coeff; - auto model = _model; - - const auto vol = particles.sliceVolume(); - auto theta = particles.sliceDilatation(); - auto m = particles.sliceWeightedVolume(); - auto force_full = KOKKOS_LAMBDA( const int i, const int j ) - { - double fx_i = 0.0; - double fy_i = 0.0; - double fz_i = 0.0; - - double xi, r, s; - double rx, ry, rz; - getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz ); - 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 ); - fx_i = coeff * rx / r; - fy_i = coeff * ry / r; - fz_i = coeff * rz / r; - - f( i, 0 ) += fx_i; - f( i, 1 ) += fy_i; - f( i, 2 ) += fz_i; - }; - - Kokkos::RangePolicy policy( 0, n_local ); - Cabana::neighbor_parallel_for( - policy, force_full, neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::ForceLPS::computeFull" ); - } - - template - double computeEnergyFull( WType& W, const PosType& x, const PosType& u, - const ParticleType& particles, - const NeighListType& neigh_list, - const int n_local, - ParallelType& neigh_op_tag ) const - { - auto theta_coeff = _model.theta_coeff; - auto s_coeff = _model.s_coeff; - auto model = _model; - - const auto vol = particles.sliceVolume(); - const auto theta = particles.sliceDilatation(); - auto m = particles.sliceWeightedVolume(); - - auto energy_full = - KOKKOS_LAMBDA( const int i, const int j, double& Phi ) - { - double xi, r, s; - getDistance( x, u, i, j, xi, r, s ); - - double num_neighbors = static_cast( - Cabana::NeighborList::numNeighbor( neigh_list, - i ) ); - - double w = ( 1.0 / num_neighbors ) * 0.5 * theta_coeff / 3.0 * - ( theta( i ) * theta( i ) ) + - 0.5 * ( s_coeff / m( i ) ) * - model.influence_function( xi ) * s * s * xi * xi * - vol( j ); - W( i ) += w; - Phi += w * vol( i ); - }; - - double strain_energy = 0.0; - - Kokkos::RangePolicy policy( 0, n_local ); - Cabana::neighbor_parallel_reduce( - policy, energy_full, neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, strain_energy, - "CabanaPD::ForceLPS::computeEnergyFull" ); - - return strain_energy; - } -}; - -template -class Force> - : public Force> -{ - protected: - using base_type = Force>; - using base_type::_half_neigh; - ForceModel _model; - - public: - using exec_space = ExecutionSpace; - - Force( const bool half_neigh, const ForceModel model ) - : base_type( half_neigh, model ) - , _model( model ) - { - } - - template - void computeWeightedVolume( ParticleType& particles, - const NeighListType& neigh_list, - const MuView& mu ) - { - auto n_local = particles.n_local; - auto x = particles.sliceReferencePosition(); - auto u = particles.sliceDisplacement(); - const auto vol = particles.sliceVolume(); - auto m = particles.sliceWeightedVolume(); - Cabana::deep_copy( m, 0.0 ); - auto model = _model; - - auto weighted_volume = KOKKOS_LAMBDA( const int i ) - { - 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 ); - - // 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 ); - } - - template - void computeDilatation( ParticleType& particles, - const NeighListType& neigh_list, - const MuView& mu ) const - { - auto n_local = particles.n_local; - const auto x = particles.sliceReferencePosition(); - auto u = particles.sliceDisplacement(); - const auto vol = particles.sliceVolume(); - auto m = particles.sliceWeightedVolume(); - auto theta = particles.sliceDilatation(); - auto model = _model; - Cabana::deep_copy( theta, 0.0 ); - - auto dilatation = KOKKOS_LAMBDA( const int i ) - { - 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 ); - - // 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, because 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 ); - } - - 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 - { - auto break_coeff = _model.bond_break_coeff; - auto theta_coeff = _model.theta_coeff; - auto s_coeff = _model.s_coeff; - auto model = _model; - - const auto vol = particles.sliceVolume(); - auto theta = particles.sliceDilatation(); - auto m = particles.sliceWeightedVolume(); - const auto nofail = particles.sliceNoFail(); - auto force_full = KOKKOS_LAMBDA( const int i ) - { - 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; - - 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; - } - } - }; - - Kokkos::RangePolicy policy( 0, n_local ); - Kokkos::parallel_for( "CabanaPD::ForceLPSDamage::computeFull", policy, - force_full ); - } - - template - double computeEnergyFull( WType& W, const PosType& x, const PosType& u, - DamageType& phi, const ParticleType& particles, - const NeighListType& neigh_list, MuView& mu, - const int n_local, ParallelType& ) const - { - auto theta_coeff = _model.theta_coeff; - auto s_coeff = _model.s_coeff; - auto model = _model; - - const auto vol = particles.sliceVolume(); - const auto theta = particles.sliceDilatation(); - auto m = particles.sliceWeightedVolume(); - - auto energy_full = KOKKOS_LAMBDA( const int i, double& Phi ) - { - std::size_t num_neighbors = - Cabana::NeighborList::numNeighbor( neigh_list, - i ); - double num_bonds = 0.0; - for ( std::size_t n = 0; n < num_neighbors; n++ ) - num_bonds += static_cast( mu( i, n ) ); - - double phi_i = 0.0; - double vol_H_i = 0.0; - for ( std::size_t n = 0; n < num_neighbors; 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 ); - - double w = - mu( i, n ) * ( ( 1.0 / num_bonds ) * 0.5 * theta_coeff / - 3.0 * ( theta( i ) * theta( i ) ) + - 0.5 * ( s_coeff / m( i ) ) * - model.influence_function( xi ) * s * s * - xi * xi * vol( j ) ); - W( i ) += w; - - phi_i += mu( i, n ) * vol( j ); - vol_H_i += vol( j ); - } - Phi += W( i ) * vol( i ); - phi( i ) = 1 - phi_i / vol_H_i; - }; - - double strain_energy = 0.0; - Kokkos::RangePolicy policy( 0, n_local ); - Kokkos::parallel_reduce( "CabanaPD::ForceLPSDamage::computeEnergyFull", - policy, energy_full, strain_energy ); - - return strain_energy; - } -}; - -template -class Force> - : public Force> -{ - protected: - using base_type = Force>; - using base_type::_half_neigh; - ForceModel _model; - - public: - using exec_space = ExecutionSpace; - - Force( const bool half_neigh, const ForceModel model ) - : base_type( half_neigh, model ) - , _model( model ) - { - } - - template - void computeForceFull( ForceType& f, const PosType& x, const PosType& u, - const ParticleType& particles, - const NeighListType& neigh_list, const int n_local, - ParallelType& neigh_op_tag ) const - { - auto theta_coeff = _model.theta_coeff; - auto s_coeff = _model.s_coeff; - auto model = _model; - - const auto vol = particles.sliceVolume(); - auto linear_theta = particles.sliceDilatation(); - // Using weighted volume from base LPS class. - auto m = particles.sliceWeightedVolume(); - - auto force_full = KOKKOS_LAMBDA( const int i, const int j ) - { - double fx_i = 0.0; - double fy_i = 0.0; - double fz_i = 0.0; - - // Get the bond distance and linearized stretch. - double xi, linear_s; - double xi_x, xi_y, xi_z; - getLinearizedDistanceComponents( x, u, i, j, xi, linear_s, xi_x, - xi_y, xi_z ); - - const double coeff = - ( theta_coeff * ( linear_theta( i ) / m( i ) + - linear_theta( j ) / m( j ) ) + - s_coeff * linear_s * ( 1.0 / m( i ) + 1.0 / m( j ) ) ) * - model.influence_function( xi ) * xi * vol( j ); - fx_i = coeff * xi_x / xi; - fy_i = coeff * xi_y / xi; - fz_i = coeff * xi_z / xi; - - f( i, 0 ) += fx_i; - f( i, 1 ) += fy_i; - f( i, 2 ) += fz_i; - }; - - Kokkos::RangePolicy policy( 0, n_local ); - Cabana::neighbor_parallel_for( - policy, force_full, neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::ForceLPS::computeFull" ); - } - - template - double computeEnergyFull( WType& W, const PosType& x, const PosType& u, - const ParticleType& particles, - const NeighListType& neigh_list, - const int n_local, - ParallelType& neigh_op_tag ) const - { - auto theta_coeff = _model.theta_coeff; - auto s_coeff = _model.s_coeff; - auto model = _model; - - const auto vol = particles.sliceVolume(); - const auto linear_theta = particles.sliceDilatation(); - auto m = particles.sliceWeightedVolume(); - - auto energy_full = - KOKKOS_LAMBDA( const int i, const int j, double& Phi ) - { - // Do we need to recompute linear_theta_i? - - double xi, linear_s; - getLinearizedDistance( x, u, i, j, xi, linear_s ); - - double num_neighbors = static_cast( - Cabana::NeighborList::numNeighbor( neigh_list, - i ) ); - - double w = ( 1.0 / num_neighbors ) * 0.5 * theta_coeff / 3.0 * - ( linear_theta( i ) * linear_theta( i ) ) + - 0.5 * ( s_coeff / m( i ) ) * - model.influence_function( xi ) * linear_s * - linear_s * xi * xi * vol( j ); - W( i ) += w; - Phi += w * vol( i ); - }; - - double strain_energy = 0.0; - - Kokkos::RangePolicy policy( 0, n_local ); - Cabana::neighbor_parallel_reduce( - policy, energy_full, neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, strain_energy, - "CabanaPD::ForceLPS::computeEnergyFull" ); - - return strain_energy; - } -}; - /****************************************************************************** - Force computation PMB + Force free functions. ******************************************************************************/ -template -class Force> -{ - protected: - bool _half_neigh; - ForceModel _model; - - public: - using exec_space = ExecutionSpace; - - Force( const bool half_neigh, const ForceModel model ) - : _half_neigh( half_neigh ) - , _model( model ) - { - } - - template - void computeWeightedVolume( ParticleType&, const NeighListType&, - const ParallelType ) - { - } - template - void computeDilatation( ParticleType&, const NeighListType&, - const ParallelType ) - { - } - - template - void computeForceFull( ForceType& f, const PosType& x, const PosType& u, - const ParticleType& particles, - const NeighListType& neigh_list, const int n_local, - ParallelType& neigh_op_tag ) const - { - auto c = _model.c; - const auto vol = particles.sliceVolume(); - - auto force_full = KOKKOS_LAMBDA( const int i, const int j ) - { - double fx_i = 0.0; - double fy_i = 0.0; - double fz_i = 0.0; - - double xi, r, s; - double rx, ry, rz; - getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz ); - const double coeff = c * s * vol( j ); - fx_i = coeff * rx / r; - fy_i = coeff * ry / r; - fz_i = coeff * rz / r; - - f( i, 0 ) += fx_i; - f( i, 1 ) += fy_i; - f( i, 2 ) += fz_i; - }; - - Kokkos::RangePolicy policy( 0, n_local ); - Cabana::neighbor_parallel_for( - policy, force_full, neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::ForcePMB::computeFull" ); - } - - template - double computeEnergyFull( WType& W, const PosType& x, const PosType& u, - const ParticleType& particles, - const NeighListType& neigh_list, - const int n_local, - ParallelType& neigh_op_tag ) const - { - auto c = _model.c; - const auto vol = particles.sliceVolume(); - - auto energy_full = - KOKKOS_LAMBDA( const int i, const int j, double& Phi ) - { - // Get the bond distance, displacement, and stretch. - double xi, r, s; - getDistance( x, u, i, j, xi, r, s ); - - // 0.25 factor is due to 1/2 from outside the integral and 1/2 from - // the integrand (pairwise potential). - double w = 0.25 * c * s * s * xi * vol( j ); - W( i ) += w; - Phi += w * vol( i ); - }; - - double strain_energy = 0.0; - Kokkos::RangePolicy policy( 0, n_local ); - Cabana::neighbor_parallel_reduce( - policy, energy_full, neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, strain_energy, - "CabanaPD::ForcePMB::computeEnergyFull" ); - - return strain_energy; - } -}; - -template -class Force> - : public Force> -{ - protected: - using base_type = Force>; - using base_type::_half_neigh; - ForceModel _model; - - public: - using exec_space = ExecutionSpace; - - Force( const bool half_neigh, const ForceModel model ) - : base_type( half_neigh, model ) - , _model( model ) - { - } - - 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 - { - 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 ) - { - 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; - - 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; - } - } - }; - - Kokkos::RangePolicy policy( 0, n_local ); - Kokkos::parallel_for( "CabanaPD::ForcePMBDamage::computeFull", policy, - force_full ); - } - - template - double computeEnergyFull( WType& W, const PosType& x, const PosType& u, - DamageType& phi, const ParticleType& particles, - const NeighListType& neigh_list, MuView& mu, - const int n_local, ParallelType& ) const - { - auto c = _model.c; - const auto vol = particles.sliceVolume(); - - auto energy_full = KOKKOS_LAMBDA( const int i, double& Phi ) - { - std::size_t num_neighbors = - Cabana::NeighborList::numNeighbor( neigh_list, - i ); - double phi_i = 0.0; - double vol_H_i = 0.0; - for ( std::size_t n = 0; n < num_neighbors; 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 ); - - // 0.25 factor is due to 1/2 from outside the integral and 1/2 - // from the integrand (pairwise potential). - double w = mu( i, n ) * 0.25 * c * s * s * xi * vol( j ); - W( i ) += w; - - phi_i += mu( i, n ) * vol( j ); - vol_H_i += vol( j ); - } - Phi += W( i ) * vol( i ); - phi( i ) = 1 - phi_i / vol_H_i; - }; - - double strain_energy = 0.0; - Kokkos::RangePolicy policy( 0, n_local ); - Kokkos::parallel_reduce( "CabanaPD::ForcePMBDamage::computeEnergyFull", - policy, energy_full, strain_energy ); - - return strain_energy; - } -}; - -template -class Force> - : public Force> -{ - protected: - using base_type = Force>; - using base_type::_half_neigh; - ForceModel _model; - - public: - using exec_space = ExecutionSpace; - - Force( const bool half_neigh, const ForceModel model ) - : base_type( half_neigh, model ) - , _model( model ) - { - } - - template - void computeForceFull( ForceType& f, const PosType& x, const PosType& u, - ParticleType& particles, - const NeighListType& neigh_list, const int n_local, - ParallelType& neigh_op_tag ) const - { - auto c = _model.c; - const auto vol = particles.sliceVolume(); - - auto force_full = KOKKOS_LAMBDA( const int i, const int j ) - { - double fx_i = 0.0; - double fy_i = 0.0; - double fz_i = 0.0; - - // Get the bond distance, displacement, and linearized stretch. - double xi, linear_s; - double xi_x, xi_y, xi_z; - getLinearizedDistanceComponents( x, u, i, j, xi, linear_s, xi_x, - xi_y, xi_z ); - - const double coeff = c * linear_s * vol( j ); - fx_i = coeff * xi_x / xi; - fy_i = coeff * xi_y / xi; - fz_i = coeff * xi_z / xi; - - f( i, 0 ) += fx_i; - f( i, 1 ) += fy_i; - f( i, 2 ) += fz_i; - }; - - Kokkos::RangePolicy policy( 0, n_local ); - Cabana::neighbor_parallel_for( - policy, force_full, neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::ForceLinearPMB::computeFull" ); - } - - template - double - computeEnergyFull( WType& W, const PosType& x, const PosType& u, - ParticleType& particles, const NeighListType& neigh_list, - const int n_local, ParallelType& neigh_op_tag ) const - { - auto c = _model.c; - const auto vol = particles.sliceVolume(); - - auto energy_full = - KOKKOS_LAMBDA( const int i, const int j, double& Phi ) - { - // Get the bond distance, displacement, and linearized stretch. - double xi, linear_s; - getLinearizedDistance( x, u, i, j, xi, linear_s ); - - // 0.25 factor is due to 1/2 from outside the integral and 1/2 from - // the integrand (pairwise potential). - double w = 0.25 * c * linear_s * linear_s * xi * vol( j ); - W( i ) += w; - Phi += w * vol( i ); - }; - - double strain_energy = 0.0; - Kokkos::RangePolicy policy( 0, n_local ); - Cabana::neighbor_parallel_reduce( - policy, energy_full, neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, strain_energy, - "CabanaPD::ForceLinearPMB::computeEnergyFull" ); - - return strain_energy; - } -}; - template void computeForce( const ForceType& force, ParticleType& particles, diff --git a/src/CabanaPD_ForceModels.hpp b/src/CabanaPD_ForceModels.hpp new file mode 100644 index 00000000..af1cf283 --- /dev/null +++ b/src/CabanaPD_ForceModels.hpp @@ -0,0 +1,33 @@ +/**************************************************************************** + * Copyright (c) 2022-2023 by Oak Ridge National Laboratory * + * All rights reserved. * + * * + * This file is part of CabanaPD. CabanaPD is distributed under a * + * BSD 3-clause license. For the licensing terms see the LICENSE file in * + * the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +#ifndef FORCE_MODELS_H +#define FORCE_MODELS_H + +namespace CabanaPD +{ +struct BaseForceModel +{ + double delta; + + BaseForceModel(){}; + BaseForceModel( const double _delta ) + : delta( _delta ){}; + + BaseForceModel( const BaseForceModel& model ) { delta = model.delta; } +}; + +template +struct ForceModel; + +} // namespace CabanaPD + +#endif diff --git a/src/force/CabanaPD_ForceModels_LPS.hpp b/src/force/CabanaPD_ForceModels_LPS.hpp new file mode 100644 index 00000000..3c0cc507 --- /dev/null +++ b/src/force/CabanaPD_ForceModels_LPS.hpp @@ -0,0 +1,146 @@ +/**************************************************************************** + * Copyright (c) 2022-2023 by Oak Ridge National Laboratory * + * All rights reserved. * + * * + * This file is part of CabanaPD. CabanaPD is distributed under a * + * BSD 3-clause license. For the licensing terms see the LICENSE file in * + * the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +#ifndef FORCE_MODELS_LPS_H +#define FORCE_MODELS_LPS_H + +#include +#include + +namespace CabanaPD +{ +template <> +struct ForceModel : public BaseForceModel +{ + using base_type = BaseForceModel; + using base_model = LPS; + using fracture_type = Elastic; + + using base_type::delta; + + int influence_type; + + double K; + double G; + double theta_coeff; + double s_coeff; + + ForceModel(){}; + ForceModel( const double _delta, const double _K, const double _G, + const int _influence = 0 ) + : base_type( _delta ) + , influence_type( _influence ) + { + set_param( _delta, _K, _G ); + } + + void set_param( const double _delta, const double _K, const double _G ) + { + delta = _delta; + K = _K; + G = _G; + + theta_coeff = 3.0 * K - 5.0 * G; + s_coeff = 15.0 * G; + } + + KOKKOS_INLINE_FUNCTION double influence_function( double xi ) const + { + if ( influence_type == 1 ) + return 1.0 / xi; + else + return 1.0; + } +}; + +template <> +struct ForceModel : public ForceModel +{ + using base_type = ForceModel; + using base_model = typename base_type::base_model; + using fracture_type = Fracture; + + using base_type::delta; + using base_type::G; + using base_type::influence_type; + using base_type::K; + using base_type::s_coeff; + using base_type::theta_coeff; + double G0; + double s0; + double bond_break_coeff; + + ForceModel() {} + ForceModel( const double _delta, const double _K, const double _G, + const double _G0, const int _influence = 0 ) + : base_type( _delta, _K, _G, _influence ) + { + set_param( _delta, _K, _G, _G0 ); + } + + void set_param( const double _delta, const double _K, const double _G, + const double _G0 ) + { + base_type::set_param( _delta, _K, _G ); + G0 = _G0; + if ( influence_type == 1 ) + { + s0 = sqrt( 5.0 * G0 / 9.0 / K / delta ); // 1/xi + } + else + { + s0 = sqrt( 8.0 * G0 / 15.0 / K / delta ); // 1 + } + bond_break_coeff = ( 1.0 + s0 ) * ( 1.0 + s0 ); + } +}; + +template <> +struct ForceModel : public ForceModel +{ + using base_type = ForceModel; + using base_model = typename base_type::base_model; + using fracture_type = typename base_type::fracture_type; + + using base_type::base_type; + + using base_type::delta; + using base_type::G; + using base_type::influence_type; + using base_type::K; + using base_type::s_coeff; + using base_type::theta_coeff; +}; + +template <> +struct ForceModel : public ForceModel +{ + using base_type = ForceModel; + using base_model = typename base_type::base_model; + using fracture_type = typename base_type::fracture_type; + + using base_type::base_type; + + using base_type::delta; + using base_type::G; + using base_type::influence_type; + using base_type::K; + using base_type::s_coeff; + using base_type::theta_coeff; + + using base_type::bond_break_coeff; + using base_type::G0; + using base_type::s0; +}; + +} // namespace CabanaPD + +#endif diff --git a/src/force/CabanaPD_ForceModels_PMB.hpp b/src/force/CabanaPD_ForceModels_PMB.hpp new file mode 100644 index 00000000..3d36740a --- /dev/null +++ b/src/force/CabanaPD_ForceModels_PMB.hpp @@ -0,0 +1,126 @@ +/**************************************************************************** + * Copyright (c) 2022-2023 by Oak Ridge National Laboratory * + * All rights reserved. * + * * + * This file is part of CabanaPD. CabanaPD is distributed under a * + * BSD 3-clause license. For the licensing terms see the LICENSE file in * + * the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +#ifndef FORCE_MODELS_PMB_H +#define FORCE_MODELS_PMB_H + +#include +#include + +namespace CabanaPD +{ +template <> +struct ForceModel : public BaseForceModel +{ + using base_type = BaseForceModel; + using base_model = PMB; + using fracture_type = Elastic; + + using base_type::delta; + + double c; + double K; + + ForceModel(){}; + ForceModel( const double delta, const double K ) + : base_type( delta ) + { + set_param( delta, K ); + } + + ForceModel( const ForceModel& model ) + : base_type( model ) + { + c = model.c; + K = model.K; + } + + void set_param( const double _delta, const double _K ) + { + delta = _delta; + K = _K; + c = 18.0 * K / ( 3.141592653589793 * delta * delta * delta * delta ); + } +}; + +template <> +struct ForceModel : public ForceModel +{ + using base_type = ForceModel; + using base_model = typename base_type::base_model; + using fracture_type = Fracture; + + using base_type::c; + using base_type::delta; + using base_type::K; + double G0; + double s0; + double bond_break_coeff; + + ForceModel() {} + ForceModel( const double delta, const double K, const double G0 ) + : base_type( delta, K ) + { + set_param( delta, K, G0 ); + } + + ForceModel( const ForceModel& model ) + : base_type( model ) + { + G0 = model.G0; + s0 = model.s0; + bond_break_coeff = model.bond_break_coeff; + } + + void set_param( const double _delta, const double _K, const double _G0 ) + { + base_type::set_param( _delta, _K ); + G0 = _G0; + s0 = sqrt( 5.0 * G0 / 9.0 / K / delta ); + bond_break_coeff = ( 1.0 + s0 ) * ( 1.0 + s0 ); + } +}; + +template <> +struct ForceModel : public ForceModel +{ + using base_type = ForceModel; + using base_model = typename base_type::base_model; + using fracture_type = typename base_type::fracture_type; + + using base_type::base_type; + + using base_type::c; + using base_type::delta; + using base_type::K; +}; + +template <> +struct ForceModel : public ForceModel +{ + using base_type = ForceModel; + using base_model = typename base_type::base_model; + using fracture_type = typename base_type::fracture_type; + + using base_type::base_type; + + using base_type::c; + using base_type::delta; + using base_type::K; + + using base_type::bond_break_coeff; + using base_type::G0; + using base_type::s0; +}; + +} // namespace CabanaPD + +#endif diff --git a/src/force/CabanaPD_Force_LPS.hpp b/src/force/CabanaPD_Force_LPS.hpp new file mode 100644 index 00000000..ffc99fa4 --- /dev/null +++ b/src/force/CabanaPD_Force_LPS.hpp @@ -0,0 +1,581 @@ +/**************************************************************************** + * Copyright (c) 2022-2023 by Oak Ridge National Laboratory * + * All rights reserved. * + * * + * This file is part of CabanaPD. CabanaPD is distributed under a * + * BSD 3-clause license. For the licensing terms see the LICENSE file in * + * the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +/**************************************************************************** + * Copyright (c) 2018-2021 by the Cabana authors * + * All rights reserved. * + * * + * This file is part of the Cabana library. Cabana is distributed under a * + * BSD 3-clause license. For the licensing terms see the LICENSE file in * + * the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +//************************************************************************ +// ExaMiniMD v. 1.0 +// Copyright (2018) National Technology & Engineering Solutions of Sandia, +// LLC (NTESS). +// +// Under the terms of Contract DE-NA-0003525 with NTESS, the U.S. Government +// retains certain rights in this software. +// +// ExaMiniMD is licensed under 3-clause BSD terms of use: Redistribution and +// use in source and binary forms, with or without modification, are +// permitted provided that the following conditions are met: +// +// 1. Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the contributors +// may be used to endorse or promote products derived from this software +// without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY NTESS "AS IS" AND ANY EXPRESS OR IMPLIED +// WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +// MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. +// IN NO EVENT SHALL NTESS OR THE CONTRIBUTORS BE LIABLE FOR ANY DIRECT, +// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +// HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, +// STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING +// IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// +//************************************************************************ + +#ifndef FORCE_LPS_H +#define FORCE_LPS_H + +#include + +#include +#include +#include +#include + +namespace CabanaPD +{ +template +class Force> +{ + protected: + bool _half_neigh; + ForceModel _model; + + public: + using exec_space = ExecutionSpace; + + Force( const bool half_neigh, const ForceModel model ) + : _half_neigh( half_neigh ) + , _model( model ) + { + } + + template + void computeWeightedVolume( ParticleType& particles, + const NeighListType& neigh_list, + const ParallelType neigh_op_tag ) + { + auto n_local = particles.n_local; + auto x = particles.sliceReferencePosition(); + auto u = particles.sliceDisplacement(); + const auto vol = particles.sliceVolume(); + auto m = particles.sliceWeightedVolume(); + Cabana::deep_copy( m, 0.0 ); + auto model = _model; + + auto weighted_volume = KOKKOS_LAMBDA( const int i, const int j ) + { + // Get the reference positions and displacements. + double xi, r, s; + getDistance( x, u, i, j, xi, r, s ); + double m_j = model.influence_function( xi ) * xi * xi * vol( j ); + m( i ) += m_j; + }; + + Kokkos::RangePolicy policy( 0, n_local ); + Cabana::neighbor_parallel_for( + policy, weighted_volume, neigh_list, Cabana::FirstNeighborsTag(), + neigh_op_tag, "CabanaPD::ForceLPS::computeWeightedVolume" ); + } + + template + void computeDilatation( ParticleType& particles, + const NeighListType& neigh_list, + const ParallelType neigh_op_tag ) const + { + auto n_local = particles.n_local; + const auto x = particles.sliceReferencePosition(); + auto u = particles.sliceDisplacement(); + const auto vol = particles.sliceVolume(); + auto m = particles.sliceWeightedVolume(); + auto theta = particles.sliceDilatation(); + auto model = _model; + Cabana::deep_copy( theta, 0.0 ); + + auto dilatation = KOKKOS_LAMBDA( const int i, const int j ) + { + // Get the bond distance, displacement, and stretch. + double xi, r, s; + getDistance( x, u, i, j, xi, r, s ); + double theta_i = + model.influence_function( xi ) * s * xi * xi * vol( j ); + theta( i ) += 3.0 * theta_i / m( i ); + }; + + Kokkos::RangePolicy policy( 0, n_local ); + Cabana::neighbor_parallel_for( + policy, dilatation, neigh_list, Cabana::FirstNeighborsTag(), + neigh_op_tag, "CabanaPD::ForceLPS::computeDilatation" ); + } + + template + void computeForceFull( ForceType& f, const PosType& x, const PosType& u, + const ParticleType& particles, + const NeighListType& neigh_list, const int n_local, + ParallelType& neigh_op_tag ) const + { + auto theta_coeff = _model.theta_coeff; + auto s_coeff = _model.s_coeff; + auto model = _model; + + const auto vol = particles.sliceVolume(); + auto theta = particles.sliceDilatation(); + auto m = particles.sliceWeightedVolume(); + auto force_full = KOKKOS_LAMBDA( const int i, const int j ) + { + double fx_i = 0.0; + double fy_i = 0.0; + double fz_i = 0.0; + + double xi, r, s; + double rx, ry, rz; + getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz ); + 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 ); + fx_i = coeff * rx / r; + fy_i = coeff * ry / r; + fz_i = coeff * rz / r; + + f( i, 0 ) += fx_i; + f( i, 1 ) += fy_i; + f( i, 2 ) += fz_i; + }; + + Kokkos::RangePolicy policy( 0, n_local ); + Cabana::neighbor_parallel_for( + policy, force_full, neigh_list, Cabana::FirstNeighborsTag(), + neigh_op_tag, "CabanaPD::ForceLPS::computeFull" ); + } + + template + double computeEnergyFull( WType& W, const PosType& x, const PosType& u, + const ParticleType& particles, + const NeighListType& neigh_list, + const int n_local, + ParallelType& neigh_op_tag ) const + { + auto theta_coeff = _model.theta_coeff; + auto s_coeff = _model.s_coeff; + auto model = _model; + + const auto vol = particles.sliceVolume(); + const auto theta = particles.sliceDilatation(); + auto m = particles.sliceWeightedVolume(); + + auto energy_full = + KOKKOS_LAMBDA( const int i, const int j, double& Phi ) + { + double xi, r, s; + getDistance( x, u, i, j, xi, r, s ); + + double num_neighbors = static_cast( + Cabana::NeighborList::numNeighbor( neigh_list, + i ) ); + + double w = ( 1.0 / num_neighbors ) * 0.5 * theta_coeff / 3.0 * + ( theta( i ) * theta( i ) ) + + 0.5 * ( s_coeff / m( i ) ) * + model.influence_function( xi ) * s * s * xi * xi * + vol( j ); + W( i ) += w; + Phi += w * vol( i ); + }; + + double strain_energy = 0.0; + + Kokkos::RangePolicy policy( 0, n_local ); + Cabana::neighbor_parallel_reduce( + policy, energy_full, neigh_list, Cabana::FirstNeighborsTag(), + neigh_op_tag, strain_energy, + "CabanaPD::ForceLPS::computeEnergyFull" ); + + return strain_energy; + } +}; + +template +class Force> + : public Force> +{ + protected: + using base_type = Force>; + using base_type::_half_neigh; + ForceModel _model; + + public: + using exec_space = ExecutionSpace; + + Force( const bool half_neigh, const ForceModel model ) + : base_type( half_neigh, model ) + , _model( model ) + { + } + + template + void computeWeightedVolume( ParticleType& particles, + const NeighListType& neigh_list, + const MuView& mu ) + { + auto n_local = particles.n_local; + auto x = particles.sliceReferencePosition(); + auto u = particles.sliceDisplacement(); + const auto vol = particles.sliceVolume(); + auto m = particles.sliceWeightedVolume(); + Cabana::deep_copy( m, 0.0 ); + auto model = _model; + + auto weighted_volume = KOKKOS_LAMBDA( const int i ) + { + 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 ); + + // 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 ); + } + + template + void computeDilatation( ParticleType& particles, + const NeighListType& neigh_list, + const MuView& mu ) const + { + auto n_local = particles.n_local; + const auto x = particles.sliceReferencePosition(); + auto u = particles.sliceDisplacement(); + const auto vol = particles.sliceVolume(); + auto m = particles.sliceWeightedVolume(); + auto theta = particles.sliceDilatation(); + auto model = _model; + Cabana::deep_copy( theta, 0.0 ); + + auto dilatation = KOKKOS_LAMBDA( const int i ) + { + 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 ); + + // 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, because 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 ); + } + + 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 + { + auto break_coeff = _model.bond_break_coeff; + auto theta_coeff = _model.theta_coeff; + auto s_coeff = _model.s_coeff; + auto model = _model; + + const auto vol = particles.sliceVolume(); + auto theta = particles.sliceDilatation(); + auto m = particles.sliceWeightedVolume(); + const auto nofail = particles.sliceNoFail(); + auto force_full = KOKKOS_LAMBDA( const int i ) + { + 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; + + 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; + } + } + }; + + Kokkos::RangePolicy policy( 0, n_local ); + Kokkos::parallel_for( "CabanaPD::ForceLPSDamage::computeFull", policy, + force_full ); + } + + template + double computeEnergyFull( WType& W, const PosType& x, const PosType& u, + DamageType& phi, const ParticleType& particles, + const NeighListType& neigh_list, MuView& mu, + const int n_local, ParallelType& ) const + { + auto theta_coeff = _model.theta_coeff; + auto s_coeff = _model.s_coeff; + auto model = _model; + + const auto vol = particles.sliceVolume(); + const auto theta = particles.sliceDilatation(); + auto m = particles.sliceWeightedVolume(); + + auto energy_full = KOKKOS_LAMBDA( const int i, double& Phi ) + { + std::size_t num_neighbors = + Cabana::NeighborList::numNeighbor( neigh_list, + i ); + double num_bonds = 0.0; + for ( std::size_t n = 0; n < num_neighbors; n++ ) + num_bonds += static_cast( mu( i, n ) ); + + double phi_i = 0.0; + double vol_H_i = 0.0; + for ( std::size_t n = 0; n < num_neighbors; 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 ); + + double w = + mu( i, n ) * ( ( 1.0 / num_bonds ) * 0.5 * theta_coeff / + 3.0 * ( theta( i ) * theta( i ) ) + + 0.5 * ( s_coeff / m( i ) ) * + model.influence_function( xi ) * s * s * + xi * xi * vol( j ) ); + W( i ) += w; + + phi_i += mu( i, n ) * vol( j ); + vol_H_i += vol( j ); + } + Phi += W( i ) * vol( i ); + phi( i ) = 1 - phi_i / vol_H_i; + }; + + double strain_energy = 0.0; + Kokkos::RangePolicy policy( 0, n_local ); + Kokkos::parallel_reduce( "CabanaPD::ForceLPSDamage::computeEnergyFull", + policy, energy_full, strain_energy ); + + return strain_energy; + } +}; + +template +class Force> + : public Force> +{ + protected: + using base_type = Force>; + using base_type::_half_neigh; + ForceModel _model; + + public: + using exec_space = ExecutionSpace; + + Force( const bool half_neigh, const ForceModel model ) + : base_type( half_neigh, model ) + , _model( model ) + { + } + + template + void computeForceFull( ForceType& f, const PosType& x, const PosType& u, + const ParticleType& particles, + const NeighListType& neigh_list, const int n_local, + ParallelType& neigh_op_tag ) const + { + auto theta_coeff = _model.theta_coeff; + auto s_coeff = _model.s_coeff; + auto model = _model; + + const auto vol = particles.sliceVolume(); + auto linear_theta = particles.sliceDilatation(); + // Using weighted volume from base LPS class. + auto m = particles.sliceWeightedVolume(); + + auto force_full = KOKKOS_LAMBDA( const int i, const int j ) + { + double fx_i = 0.0; + double fy_i = 0.0; + double fz_i = 0.0; + + // Get the bond distance and linearized stretch. + double xi, linear_s; + double xi_x, xi_y, xi_z; + getLinearizedDistanceComponents( x, u, i, j, xi, linear_s, xi_x, + xi_y, xi_z ); + + const double coeff = + ( theta_coeff * ( linear_theta( i ) / m( i ) + + linear_theta( j ) / m( j ) ) + + s_coeff * linear_s * ( 1.0 / m( i ) + 1.0 / m( j ) ) ) * + model.influence_function( xi ) * xi * vol( j ); + fx_i = coeff * xi_x / xi; + fy_i = coeff * xi_y / xi; + fz_i = coeff * xi_z / xi; + + f( i, 0 ) += fx_i; + f( i, 1 ) += fy_i; + f( i, 2 ) += fz_i; + }; + + Kokkos::RangePolicy policy( 0, n_local ); + Cabana::neighbor_parallel_for( + policy, force_full, neigh_list, Cabana::FirstNeighborsTag(), + neigh_op_tag, "CabanaPD::ForceLPS::computeFull" ); + } + + template + double computeEnergyFull( WType& W, const PosType& x, const PosType& u, + const ParticleType& particles, + const NeighListType& neigh_list, + const int n_local, + ParallelType& neigh_op_tag ) const + { + auto theta_coeff = _model.theta_coeff; + auto s_coeff = _model.s_coeff; + auto model = _model; + + const auto vol = particles.sliceVolume(); + const auto linear_theta = particles.sliceDilatation(); + auto m = particles.sliceWeightedVolume(); + + auto energy_full = + KOKKOS_LAMBDA( const int i, const int j, double& Phi ) + { + // Do we need to recompute linear_theta_i? + + double xi, linear_s; + getLinearizedDistance( x, u, i, j, xi, linear_s ); + + double num_neighbors = static_cast( + Cabana::NeighborList::numNeighbor( neigh_list, + i ) ); + + double w = ( 1.0 / num_neighbors ) * 0.5 * theta_coeff / 3.0 * + ( linear_theta( i ) * linear_theta( i ) ) + + 0.5 * ( s_coeff / m( i ) ) * + model.influence_function( xi ) * linear_s * + linear_s * xi * xi * vol( j ); + W( i ) += w; + Phi += w * vol( i ); + }; + + double strain_energy = 0.0; + + Kokkos::RangePolicy policy( 0, n_local ); + Cabana::neighbor_parallel_reduce( + policy, energy_full, neigh_list, Cabana::FirstNeighborsTag(), + neigh_op_tag, strain_energy, + "CabanaPD::ForceLPS::computeEnergyFull" ); + + return strain_energy; + } +}; + +} // namespace CabanaPD + +#endif diff --git a/src/force/CabanaPD_Force_PMB.hpp b/src/force/CabanaPD_Force_PMB.hpp new file mode 100644 index 00000000..2b1ba5dc --- /dev/null +++ b/src/force/CabanaPD_Force_PMB.hpp @@ -0,0 +1,387 @@ +/**************************************************************************** + * Copyright (c) 2022-2023 by Oak Ridge National Laboratory * + * All rights reserved. * + * * + * This file is part of CabanaPD. CabanaPD is distributed under a * + * BSD 3-clause license. For the licensing terms see the LICENSE file in * + * the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +/**************************************************************************** + * Copyright (c) 2018-2021 by the Cabana authors * + * All rights reserved. * + * * + * This file is part of the Cabana library. Cabana is distributed under a * + * BSD 3-clause license. For the licensing terms see the LICENSE file in * + * the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +//************************************************************************ +// ExaMiniMD v. 1.0 +// Copyright (2018) National Technology & Engineering Solutions of Sandia, +// LLC (NTESS). +// +// Under the terms of Contract DE-NA-0003525 with NTESS, the U.S. Government +// retains certain rights in this software. +// +// ExaMiniMD is licensed under 3-clause BSD terms of use: Redistribution and +// use in source and binary forms, with or without modification, are +// permitted provided that the following conditions are met: +// +// 1. Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the contributors +// may be used to endorse or promote products derived from this software +// without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY NTESS "AS IS" AND ANY EXPRESS OR IMPLIED +// WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +// MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. +// IN NO EVENT SHALL NTESS OR THE CONTRIBUTORS BE LIABLE FOR ANY DIRECT, +// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +// HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, +// STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING +// IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// +//************************************************************************ + +#ifndef FORCE_PMB_H +#define FORCE_PMB_H + +#include + +#include +#include +#include +#include + +namespace CabanaPD +{ +template +class Force> +{ + protected: + bool _half_neigh; + ForceModel _model; + + public: + using exec_space = ExecutionSpace; + + Force( const bool half_neigh, const ForceModel model ) + : _half_neigh( half_neigh ) + , _model( model ) + { + } + + template + void computeWeightedVolume( ParticleType&, const NeighListType&, + const ParallelType ) + { + } + template + void computeDilatation( ParticleType&, const NeighListType&, + const ParallelType ) + { + } + + template + void computeForceFull( ForceType& f, const PosType& x, const PosType& u, + const ParticleType& particles, + const NeighListType& neigh_list, const int n_local, + ParallelType& neigh_op_tag ) const + { + auto c = _model.c; + const auto vol = particles.sliceVolume(); + + auto force_full = KOKKOS_LAMBDA( const int i, const int j ) + { + double fx_i = 0.0; + double fy_i = 0.0; + double fz_i = 0.0; + + double xi, r, s; + double rx, ry, rz; + getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz ); + const double coeff = c * s * vol( j ); + fx_i = coeff * rx / r; + fy_i = coeff * ry / r; + fz_i = coeff * rz / r; + + f( i, 0 ) += fx_i; + f( i, 1 ) += fy_i; + f( i, 2 ) += fz_i; + }; + + Kokkos::RangePolicy policy( 0, n_local ); + Cabana::neighbor_parallel_for( + policy, force_full, neigh_list, Cabana::FirstNeighborsTag(), + neigh_op_tag, "CabanaPD::ForcePMB::computeFull" ); + } + + template + double computeEnergyFull( WType& W, const PosType& x, const PosType& u, + const ParticleType& particles, + const NeighListType& neigh_list, + const int n_local, + ParallelType& neigh_op_tag ) const + { + auto c = _model.c; + const auto vol = particles.sliceVolume(); + + auto energy_full = + KOKKOS_LAMBDA( const int i, const int j, double& Phi ) + { + // Get the bond distance, displacement, and stretch. + double xi, r, s; + getDistance( x, u, i, j, xi, r, s ); + + // 0.25 factor is due to 1/2 from outside the integral and 1/2 from + // the integrand (pairwise potential). + double w = 0.25 * c * s * s * xi * vol( j ); + W( i ) += w; + Phi += w * vol( i ); + }; + + double strain_energy = 0.0; + Kokkos::RangePolicy policy( 0, n_local ); + Cabana::neighbor_parallel_reduce( + policy, energy_full, neigh_list, Cabana::FirstNeighborsTag(), + neigh_op_tag, strain_energy, + "CabanaPD::ForcePMB::computeEnergyFull" ); + + return strain_energy; + } +}; + +template +class Force> + : public Force> +{ + protected: + using base_type = Force>; + using base_type::_half_neigh; + ForceModel _model; + + public: + using exec_space = ExecutionSpace; + + Force( const bool half_neigh, const ForceModel model ) + : base_type( half_neigh, model ) + , _model( model ) + { + } + + 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 + { + 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 ) + { + 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; + + 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; + } + } + }; + + Kokkos::RangePolicy policy( 0, n_local ); + Kokkos::parallel_for( "CabanaPD::ForcePMBDamage::computeFull", policy, + force_full ); + } + + template + double computeEnergyFull( WType& W, const PosType& x, const PosType& u, + DamageType& phi, const ParticleType& particles, + const NeighListType& neigh_list, MuView& mu, + const int n_local, ParallelType& ) const + { + auto c = _model.c; + const auto vol = particles.sliceVolume(); + + auto energy_full = KOKKOS_LAMBDA( const int i, double& Phi ) + { + std::size_t num_neighbors = + Cabana::NeighborList::numNeighbor( neigh_list, + i ); + double phi_i = 0.0; + double vol_H_i = 0.0; + for ( std::size_t n = 0; n < num_neighbors; 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 ); + + // 0.25 factor is due to 1/2 from outside the integral and 1/2 + // from the integrand (pairwise potential). + double w = mu( i, n ) * 0.25 * c * s * s * xi * vol( j ); + W( i ) += w; + + phi_i += mu( i, n ) * vol( j ); + vol_H_i += vol( j ); + } + Phi += W( i ) * vol( i ); + phi( i ) = 1 - phi_i / vol_H_i; + }; + + double strain_energy = 0.0; + Kokkos::RangePolicy policy( 0, n_local ); + Kokkos::parallel_reduce( "CabanaPD::ForcePMBDamage::computeEnergyFull", + policy, energy_full, strain_energy ); + + return strain_energy; + } +}; + +template +class Force> + : public Force> +{ + protected: + using base_type = Force>; + using base_type::_half_neigh; + ForceModel _model; + + public: + using exec_space = ExecutionSpace; + + Force( const bool half_neigh, const ForceModel model ) + : base_type( half_neigh, model ) + , _model( model ) + { + } + + template + void computeForceFull( ForceType& f, const PosType& x, const PosType& u, + ParticleType& particles, + const NeighListType& neigh_list, const int n_local, + ParallelType& neigh_op_tag ) const + { + auto c = _model.c; + const auto vol = particles.sliceVolume(); + + auto force_full = KOKKOS_LAMBDA( const int i, const int j ) + { + double fx_i = 0.0; + double fy_i = 0.0; + double fz_i = 0.0; + + // Get the bond distance, displacement, and linearized stretch. + double xi, linear_s; + double xi_x, xi_y, xi_z; + getLinearizedDistanceComponents( x, u, i, j, xi, linear_s, xi_x, + xi_y, xi_z ); + + const double coeff = c * linear_s * vol( j ); + fx_i = coeff * xi_x / xi; + fy_i = coeff * xi_y / xi; + fz_i = coeff * xi_z / xi; + + f( i, 0 ) += fx_i; + f( i, 1 ) += fy_i; + f( i, 2 ) += fz_i; + }; + + Kokkos::RangePolicy policy( 0, n_local ); + Cabana::neighbor_parallel_for( + policy, force_full, neigh_list, Cabana::FirstNeighborsTag(), + neigh_op_tag, "CabanaPD::ForceLinearPMB::computeFull" ); + } + + template + double + computeEnergyFull( WType& W, const PosType& x, const PosType& u, + ParticleType& particles, const NeighListType& neigh_list, + const int n_local, ParallelType& neigh_op_tag ) const + { + auto c = _model.c; + const auto vol = particles.sliceVolume(); + + auto energy_full = + KOKKOS_LAMBDA( const int i, const int j, double& Phi ) + { + // Get the bond distance, displacement, and linearized stretch. + double xi, linear_s; + getLinearizedDistance( x, u, i, j, xi, linear_s ); + + // 0.25 factor is due to 1/2 from outside the integral and 1/2 from + // the integrand (pairwise potential). + double w = 0.25 * c * linear_s * linear_s * xi * vol( j ); + W( i ) += w; + Phi += w * vol( i ); + }; + + double strain_energy = 0.0; + Kokkos::RangePolicy policy( 0, n_local ); + Cabana::neighbor_parallel_reduce( + policy, energy_full, neigh_list, Cabana::FirstNeighborsTag(), + neigh_op_tag, strain_energy, + "CabanaPD::ForceLinearPMB::computeEnergyFull" ); + + return strain_energy; + } +}; + +} // namespace CabanaPD + +#endif diff --git a/unit_test/tstForce.hpp b/unit_test/tstForce.hpp index 4aa16454..b3fc96a6 100644 --- a/unit_test/tstForce.hpp +++ b/unit_test/tstForce.hpp @@ -28,8 +28,13 @@ #include #include +#include #include #include +#include +#include +#include +#include namespace Test {