diff --git a/src/CabanaPD.hpp b/src/CabanaPD.hpp index 7da7ad22..35901f00 100644 --- a/src/CabanaPD.hpp +++ b/src/CabanaPD.hpp @@ -29,8 +29,10 @@ #include #include +#include #include #include +#include #include #include diff --git a/src/CabanaPD_ForceModels.hpp b/src/CabanaPD_ForceModels.hpp index 6685de1c..5aa8d18f 100644 --- a/src/CabanaPD_ForceModels.hpp +++ b/src/CabanaPD_ForceModels.hpp @@ -31,6 +31,29 @@ struct BaseForceModel void thermalStretch( double&, const int, const int ) const {} }; +struct BaseInfluenceForceModel : public BaseForceModel +{ + using base_type = BaseForceModel; + + using base_type::delta; + int influence_type; + + BaseInfluenceForceModel( const double _delta, const int _influence = 0 ) + : base_type( _delta ) + , influence_type( _influence ) + { + } + + KOKKOS_INLINE_FUNCTION + double influence_function( const double xi ) const + { + if ( influence_type == 1 ) + return 1.0 / xi; + else + return 1.0; + } +}; + template struct BaseTemperatureModel { diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index c9e9816d..459ddfec 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -652,6 +652,155 @@ class Particles using base_type::_timer; }; +template +class Particles + : public Particles +{ + public: + using self_type = Particles; + using base_type = + Particles; + using memory_space = typename base_type::memory_space; + using base_type::dim; + + // Per particle. + using base_type::n_ghost; + using base_type::n_global; + using base_type::n_local; + using base_type::size; + + // Split? + using tensor_type = Cabana::MemberTypes; + using aosoa_tensor_type = Cabana::AoSoA; + + // Per type. + using base_type::n_types; + + // Simulation total domain. + using base_type::global_mesh_ext; + + // Simulation sub domain (single MPI rank). + using base_type::ghost_mesh_hi; + using base_type::ghost_mesh_lo; + using base_type::local_mesh_ext; + using base_type::local_mesh_hi; + using base_type::local_mesh_lo; + + using base_type::dx; + using base_type::local_grid; + + using base_type::halo_width; + + // Default constructor. + Particles() + : base_type() + { + _init_timer.start(); + _aosoa_F = aosoa_tensor_type( "Particle Deformation Gradient", 0 ); + _aosoa_Kinv = aosoa_tensor_type( "Particle Shape Tensor", 0 ); + _init_timer.stop(); + } + + // Constructor which initializes particles on regular grid. + template + Particles( Args&&... args ) + : base_type( std::forward( args )... ) + { + _init_timer.start(); + _aosoa_F = + aosoa_tensor_type( "Particle Deformation Gradient", n_local ); + _aosoa_Kinv = aosoa_tensor_type( "Particle Shape Tensor", n_local ); + init_corr(); + _init_timer.stop(); + } + + template + void createParticles( Args&&... args ) + { + // Forward arguments to standard or custom particle creation. + base_type::createParticles( std::forward( args )... ); + _init_timer.start(); + _aosoa_F.resize( n_local ); + _aosoa_Kinv.resize( n_local ); + _init_timer.stop(); + } + + auto sliceDefGradient() + { + return Cabana::slice<0>( _aosoa_F, "deformation_gradient" ); + } + auto sliceDefGradient() const { return sliceDefGradient(); } + auto sliceShapeTensor() + { + return Cabana::slice<0>( _aosoa_Kinv, "shape_tensor" ); + } + auto sliceShapeTensor() const { return sliceShapeTensor(); } + + void resize( int new_local, int new_ghost ) + { + base_type::resize( new_local, new_ghost ); + _timer.start(); + _aosoa_F.resize( new_local + new_ghost ); + _aosoa_Kinv.resize( new_local + new_ghost ); + _timer.stop(); + } + + void output( [[maybe_unused]] const int output_step, + [[maybe_unused]] const double output_time, + [[maybe_unused]] const bool use_reference = true ) + { + _output_timer.start(); + +#ifdef Cabana_ENABLE_HDF5 + Cabana::Experimental::HDF5ParticleOutput::writeTimeStep( + h5_config, "particles", MPI_COMM_WORLD, output_step, output_time, + n_local, base_type::getPosition( use_reference ), + base_type::sliceStrainEnergy(), base_type::sliceForce(), + base_type::sliceDisplacement(), base_type::sliceVelocity(), + base_type::sliceDamage(), sliceDefGradient() ); +#else +#ifdef Cabana_ENABLE_SILO + Cabana::Grid::Experimental::SiloParticleOutput:: + writePartialRangeTimeStep( + "particles", local_grid->globalGrid(), output_step, output_time, + 0, n_local, base_type::getPosition( use_reference ), + base_type::sliceStrainEnergy(), base_type::sliceForce(), + base_type::sliceDisplacement(), base_type::sliceVelocity(), + base_type::sliceDamage(), sliceDefGradient() ); +#else + log( std::cout, "No particle output enabled." ); +#endif +#endif + + _output_timer.stop(); + } + + friend class Comm; + friend class Comm; + friend class Comm; + + protected: + void init_corr() + { + auto F = sliceDefGradient(); + Cabana::deep_copy( F, 0.0 ); + auto Kinv = sliceShapeTensor(); + Cabana::deep_copy( Kinv, 0.0 ); + } + + aosoa_tensor_type _aosoa_F; + aosoa_tensor_type _aosoa_Kinv; + +#ifdef Cabana_ENABLE_HDF5 + using base_type::h5_config; +#endif + + using base_type::_init_timer; + using base_type::_output_timer; + using base_type::_timer; +}; + template class Particles : public Particles diff --git a/src/CabanaPD_Types.hpp b/src/CabanaPD_Types.hpp index 62076735..ea4f3cee 100644 --- a/src/CabanaPD_Types.hpp +++ b/src/CabanaPD_Types.hpp @@ -39,6 +39,9 @@ struct LPS struct LinearLPS { }; +struct Correspondence +{ +}; } // namespace CabanaPD #endif diff --git a/src/force/CabanaPD_ForceModels_Correspondence.hpp b/src/force/CabanaPD_ForceModels_Correspondence.hpp new file mode 100644 index 00000000..a93e4542 --- /dev/null +++ b/src/force/CabanaPD_ForceModels_Correspondence.hpp @@ -0,0 +1,91 @@ +/**************************************************************************** + * Copyright (c) 2022 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_CORR_H +#define FORCE_MODELS_CORR_H + +#include +#include + +namespace CabanaPD +{ +template <> +struct ForceModel : public BaseInfluenceForceModel +{ + using base_type = BaseInfluenceForceModel; + using base_model = Correspondence; + using fracture_type = Elastic; + using thermal_type = TemperatureIndependent; + + using base_type::delta; + using base_type::influence_type; + + double K; + double G; + double theta_coeff; + double s_coeff; + + ForceModel( const double _delta, const double _K, const double _G, + const int _influence = 0 ) + : base_type( _delta, _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; + } + + template + KOKKOS_INLINE_FUNCTION auto getStress( DefGradType F, const int i ) + { + double epsilon[3][3]; + + for ( std::size_t d = 0; d < 3; d++ ) + { + epsilon[d][d] = F( i, d, d ) - 1.0; + int d2 = getComponent( d ); + epsilon[d][d2] = 0.5 * ( F( i, d, d2 ) + F( i, d2, d ) ); + epsilon[d2][d] = epsilon[d][d2]; + } + + double sigma[3][3]; + double diag = ( K - 2.0 / 3.0 * G ) * + ( epsilon[0][0] + epsilon[1][1] + epsilon[2][2] ); + for ( std::size_t d = 0; d < 3; d++ ) + { + sigma[d][d] = diag + 2.0 * G * epsilon[d][d]; + int d2 = getComponent( d ); + sigma[d][d2] = 2.0 * G * epsilon[d][d2]; + sigma[d2][d] = sigma[d][d2]; + } + return sigma; + } + + auto getComponent( const int d ) + { + int d2; + if ( d < 2 ) + return d + 1; + else + return 0; + } +}; + +} // namespace CabanaPD + +#endif diff --git a/src/force/CabanaPD_ForceModels_LPS.hpp b/src/force/CabanaPD_ForceModels_LPS.hpp index 8f2999fc..44e56cee 100644 --- a/src/force/CabanaPD_ForceModels_LPS.hpp +++ b/src/force/CabanaPD_ForceModels_LPS.hpp @@ -18,27 +18,24 @@ namespace CabanaPD { template <> -struct ForceModel : public BaseForceModel +struct ForceModel : public BaseInfluenceForceModel { - using base_type = BaseForceModel; + using base_type = BaseInfluenceForceModel; using base_model = LPS; using fracture_type = Elastic; using thermal_type = TemperatureIndependent; using base_type::delta; - - int influence_type; + using base_type::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 ) + : base_type( _delta, _influence ) { set_param( _delta, _K, _G ); } @@ -52,14 +49,6 @@ struct ForceModel : public BaseForceModel 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 <> @@ -80,7 +69,6 @@ struct ForceModel : public ForceModel 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 ) diff --git a/src/force/CabanaPD_Force_Correspondence.hpp b/src/force/CabanaPD_Force_Correspondence.hpp new file mode 100644 index 00000000..ec718441 --- /dev/null +++ b/src/force/CabanaPD_Force_Correspondence.hpp @@ -0,0 +1,223 @@ +/**************************************************************************** + * Copyright (c) 2022 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_CORR_H +#define FORCE_CORR_H + +#include + +#include +#include +#include +#include + +namespace CabanaPD +{ +template +class Force> + : public Force +{ + protected: + using base_type = Force; + using base_type::_half_neigh; + ForceModel _model; + + using base_type::_energy_timer; + using base_type::_timer; + + public: + using exec_space = ExecutionSpace; + + Force( const bool half_neigh, + const ForceModel model ) + : base_type( half_neigh ) + , _model( model ) + { + } + + template + void computeShapeTensor( ParticleType& particles, + const NeighListType& neigh_list, + const ParallelType neigh_op_tag ) + { + _timer.start(); + + auto n_local = particles.n_local; + auto x = particles.sliceReferencePosition(); + auto u = particles.sliceDisplacement(); + const auto vol = particles.sliceVolume(); + auto K = particles.sliceShapeTensor(); + Cabana::deep_copy( K, 0.0 ); + auto model = _model; + + auto shape = 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 K_j = model.influence_function( xi ) * xi * xi * vol( j ); + K( i ) += K_j; + }; + Kokkos::RangePolicy policy( 0, n_local ); + Cabana::neighbor_parallel_for( + policy, shape, neigh_list, Cabana::FirstNeighborsTag(), + neigh_op_tag, "CabanaPD::ForceCorrespondence::computeShapeTensor" ); + Kokkos::fence(); + + auto invert = KOKKOS_LAMBDA( const int i ) + { + auto k = K( i ); + auto det = + ( k[0][0] * k[1][1] * k[2][2] - k[0][0] * k[1][2] * k[1][2] ) - + ( k[0][1] * k[0][1] * k[2][2] - k[0][1] * k[0][2] * k[1][2] ) + + ( k[0][2] * k[0][1] * k[1][2] - k[0][2] * k[0][2] * k[1][1] ); + + K( i, 0, 0 ) = ( -k[1][2] * k[1][2] + k[1][1] * k[2] k[2] ) / det; + K( i, 1, 1 ) = ( -k[0][2] * k[0][2] + k[0][0] * k[2] k[2] ) / det; + K( i, 2, 2 ) = ( -k[0][1] * k[0][1] + k[0][0] * k[1] k[1] ) / det; + }; + Kokkos::parallel_for( + "CabanaPD::ForceCorrespondence::invertShapeTensor", policy, + invert ); + + _timer.stop(); + } + + template + void computeDeformationGradient( ParticleType& particles, + const NeighListType& neigh_list, + const ParallelType neigh_op_tag ) + { + _timer.start(); + + 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::ForceCorrespondence::computeDilatation" ); + + _timer.stop(); + } + + template + void computeForceFull( ForceType& f, const PosType&, const PosType& u, + const ParticleType& particles, + const NeighListType& neigh_list, const int n_local, + ParallelType& neigh_op_tag ) + { + _timer.start(); + + auto theta_coeff = _model.theta_coeff; + auto s_coeff = _model.s_coeff; + auto model = _model; + + const auto vol = particles.sliceVolume(); + const auto Kinv = particles.sliceShapeTensor(); + const auto y = particles.sliceCurrentPosition(); + 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( y, u, i, j, xi, r, s, rx, ry, rz ); + const double coeff = model.influence_function( xi ) * + model.getStress( F, i ) * Kinv * 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::ForceCorrespondence::computeFull" ); + + _timer.stop(); + } + + 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 ) + { + _energy_timer.start(); + + 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::ForceCorrespondence::computeEnergyFull" ); + + _energy_timer.stop(); + return strain_energy; + } +}; + +} // namespace CabanaPD + +#endif