From c65cb8710a67d85c10da556649a52803929b7867 Mon Sep 17 00:00:00 2001 From: Pablo Seleson Date: Fri, 21 Jul 2023 16:55:32 -0400 Subject: [PATCH] Initial contact forces Co-authored-by: Pablo Seleson --- examples/mechanics/CMakeLists.txt | 5 +- examples/mechanics/disk_impact.cpp | 101 ++++++++++++ examples/mechanics/inputs/disk_impact.json | 16 ++ src/CabanaPD.hpp | 1 + src/CabanaPD_Contact.hpp | 181 +++++++++++++++++++++ src/CabanaPD_Particles.hpp | 7 +- src/CabanaPD_Solver.hpp | 136 ++++++++++++++++ 7 files changed, 443 insertions(+), 4 deletions(-) create mode 100644 examples/mechanics/disk_impact.cpp create mode 100644 examples/mechanics/inputs/disk_impact.json create mode 100644 src/CabanaPD_Contact.hpp diff --git a/examples/mechanics/CMakeLists.txt b/examples/mechanics/CMakeLists.txt index a10b05f1..19148fa3 100644 --- a/examples/mechanics/CMakeLists.txt +++ b/examples/mechanics/CMakeLists.txt @@ -10,4 +10,7 @@ target_link_libraries(CrackBranching LINK_PUBLIC CabanaPD) add_executable(FragmentingCylinder fragmenting_cylinder.cpp) target_link_libraries(FragmentingCylinder LINK_PUBLIC CabanaPD) -install(TARGETS ElasticWave KalthoffWinkler CrackBranching FragmentingCylinder DESTINATION ${CMAKE_INSTALL_BINDIR}) +add_executable(DiskImpact disk_impact.cpp) +target_link_libraries(DiskImpact LINK_PUBLIC CabanaPD) + +install(TARGETS ElasticWave KalthoffWinkler CrackBranching FragmentingCylinder DiskImpact DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/examples/mechanics/disk_impact.cpp b/examples/mechanics/disk_impact.cpp new file mode 100644 index 00000000..924165e3 --- /dev/null +++ b/examples/mechanics/disk_impact.cpp @@ -0,0 +1,101 @@ +#include +#include +#include + +#include "mpi.h" + +#include + +#include + +// Simulate a sphere impacting a thin cylindrical plate. +void diskImpactExample( const std::string filename ) +{ + using exec_space = Kokkos::DefaultExecutionSpace; + using memory_space = typename exec_space::memory_space; + + CabanaPD::Inputs inputs( filename ); + + double K = inputs["bulk_modulus"]; + double G0 = inputs["fracture_energy"]; + double delta = inputs["horizon"]; + delta += 1e-10; + // Choose force model type. + using model_type = CabanaPD::ForceModel; + model_type force_model( delta, K, G0 ); + + // Create particles from mesh. + auto particles = CabanaPD::createParticles( + exec_space(), inputs ); + + std::array low_corner = inputs["low_corner"]; + std::array high_corner = inputs["high_corner"]; + double global_mesh_ext = high_corner[0] - low_corner[0]; + auto x = particles->sliceReferencePosition(); + auto create_func = KOKKOS_LAMBDA( const int, const double px[3] ) + { + auto width = global_mesh_ext / 2.0; + auto r2 = px[0] * px[0] + px[1] * px[1]; + if ( r2 > width * width ) + return false; + return true; + }; + particles->createParticles( exec_space{}, create_func ); + + // Define particle initialization. + auto v = particles->sliceVelocity(); + auto f = particles->sliceForce(); + auto rho = particles->sliceDensity(); + + double rho0 = inputs["density"]; + auto init_functor = KOKKOS_LAMBDA( const int pid ) + { + rho( pid ) = rho0; + for ( int d = 0; d < 3; d++ ) + v( pid, d ) = 0.0; + }; + particles->updateParticles( exec_space{}, init_functor ); + + double dx = particles->dx[0]; + double r_c = dx * 2.0; + CabanaPD::NormalRepulsionModel contact_model( delta, r_c, K ); + + // FIXME: use createSolver to switch backend at runtime. + auto cabana_pd = CabanaPD::createSolverContact( + inputs, particles, force_model, contact_model ); + + double impact_r = inputs["impactor_radius"]; + double impact_v = inputs["impactor_velocity"]; + double impact_x = 0.0; + double impact_y = 0.0; + auto body_func = KOKKOS_LAMBDA( const int p, const double t ) + { + auto z = t * impact_v + impact_r; + double r = sqrt( ( x( p, 0 ) - impact_x ) * ( x( p, 0 ) - impact_x ) + + ( x( p, 1 ) - impact_y ) * ( x( p, 1 ) - impact_y ) + + ( x( p, 2 ) - z ) * ( x( p, 2 ) - z ) ); + if ( r < impact_r ) + { + double fmag = -1.0e17 * ( r - impact_r ) * ( r - impact_r ); + f( p, 0 ) += fmag * x( p, 0 ) / r; + f( p, 1 ) += fmag * x( p, 1 ) / r; + f( p, 2 ) += fmag * ( x( p, 2 ) - z ) / r; + } + }; + auto body = CabanaPD::createBodyTerm( body_func, true ); + + cabana_pd->init( body ); + cabana_pd->run( body ); +} + +int main( int argc, char* argv[] ) +{ + MPI_Init( &argc, &argv ); + + Kokkos::initialize( argc, argv ); + + diskImpactExample( argv[1] ); + + Kokkos::finalize(); + MPI_Finalize(); +} diff --git a/examples/mechanics/inputs/disk_impact.json b/examples/mechanics/inputs/disk_impact.json new file mode 100644 index 00000000..6a273f0b --- /dev/null +++ b/examples/mechanics/inputs/disk_impact.json @@ -0,0 +1,16 @@ +{ + "num_cells" : {"value": [149, 149, 5]}, + "low_corner" : {"value": [-0.037, -0.037, 2.5e-3], "unit": "m"}, + "high_corner" : {"value": [0.037, 0.037, 0.0], "unit": "m"}, + "density" : {"value": 2200.0, "unit": "kg/m^3"}, + "bulk_modulus" : {"value": 14.9e+9, "unit": "Pa"}, + "fracture_energy" : {"value": 10.0575, "unit": "J/m^2"}, + "horizon" : {"value": 0.0015, "unit": "m"}, + "final_time" : {"value": 2.0e-4, "unit": "s"}, + "timestep" : {"value": 1.0e-7, "unit": "s"}, + "timestep_safety_factor" : {"value": 0.85}, + "output_frequency" : {"value": 50}, + "output_reference" : {"value": false}, + "impactor_radius" : {"value": 5e-3, "unit": "m"}, + "impactor_velocity" : {"value": -100.0, "unit": "m/s"} +} diff --git a/src/CabanaPD.hpp b/src/CabanaPD.hpp index 7da7ad22..380ab056 100644 --- a/src/CabanaPD.hpp +++ b/src/CabanaPD.hpp @@ -16,6 +16,7 @@ #include #include #include +#include #include #include #include diff --git a/src/CabanaPD_Contact.hpp b/src/CabanaPD_Contact.hpp new file mode 100644 index 00000000..02f4e896 --- /dev/null +++ b/src/CabanaPD_Contact.hpp @@ -0,0 +1,181 @@ +/**************************************************************************** + * 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 CONTACT_H +#define CONTACT_H + +#include + +#include +#include +#include + +namespace CabanaPD +{ +/****************************************************************************** + Contact model +******************************************************************************/ +struct ContactModel +{ + double delta; + double Rc; + + ContactModel(){}; + // PD horizon + // Contact radius + ContactModel( const double _delta, const double _Rc ) + : delta( _delta ) + , Rc( _Rc ){}; +}; + +/* Normal repulsion */ + +struct NormalRepulsionModel : public ContactModel +{ + using ContactModel::delta; + using ContactModel::Rc; + + double c; + double K; + + NormalRepulsionModel(){}; + NormalRepulsionModel( const double delta, const double Rc, const double _K ) + : ContactModel( delta, Rc ) + , K( _K ) + { + set_param( delta, Rc, K ); + } + + void set_param( const double _delta, const double _Rc, const double _K ) + { + delta = _delta; + Rc = _Rc; + K = _K; + // This could inherit from PMB (same c) + c = 18.0 * K / ( 3.141592653589793 * delta * delta * delta * delta ); + } +}; + +template +class Contact; + +/****************************************************************************** + Normal repulsion computation +******************************************************************************/ +template +class Contact +{ + public: + using neighbor_type = + Cabana::VerletList; + + template + Contact( const bool half_neigh, const NormalRepulsionModel model, + ParticleType& particles ) + : _half_neigh( half_neigh ) + , _model( model ) + { + // Create contact neighbor list + for ( std::size_t d = 0; d < 3; d++ ) + { + mesh_min[d] = particles.ghost_mesh_lo[d]; + mesh_max[d] = particles.ghost_mesh_hi[d]; + } + const auto y = particles.sliceCurrentPosition(); + _neighbors = std::make_shared( + y, 0, particles.n_local, model.Rc, 1.0, mesh_min, mesh_max ); + } + + template + void computeContactFull( ForceType& fc, const ParticleType& particles, + const int n_local, + ParallelType& neigh_op_tag ) const + { + auto delta = _model.delta; + auto Rc = _model.Rc; + auto c = _model.c; + const auto vol = particles.sliceVolume(); + const auto x = particles.sliceReferencePosition(); + const auto u = particles.sliceDisplacement(); + const auto y = particles.sliceCurrentPosition(); + + _neighbors->build( y, 0, particles.n_local, Rc, 1.0, mesh_min, + mesh_max ); + + auto contact_full = KOKKOS_LAMBDA( const int i, const int j ) + { + double fcx_i = 0.0; + double fcy_i = 0.0; + double fcz_i = 0.0; + + double xi, r, s; + double rx, ry, rz; + getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz ); + + // Contact "stretch" + const double sc = ( r - Rc ) / delta; + + // Normal repulsion uses a 15 factor compared to the PMB force + const double coeff = 15 * c * sc * vol( j ); + fcx_i = coeff * rx / r; + fcy_i = coeff * ry / r; + fcz_i = coeff * rz / r; + + fc( i, 0 ) += fcx_i; + fc( i, 1 ) += fcy_i; + fc( i, 2 ) += fcz_i; + }; + + // FIXME: using default space for now. + using exec_space = typename MemorySpace::execution_space; + Kokkos::RangePolicy policy( 0, n_local ); + Cabana::neighbor_parallel_for( + policy, contact_full, *_neighbors, Cabana::FirstNeighborsTag(), + neigh_op_tag, "CabanaPD::Contact::compute_full" ); + } + + protected: + bool _half_neigh; + NormalRepulsionModel _model; + std::shared_ptr _neighbors; + + double mesh_max[3]; + double mesh_min[3]; +}; + +template +void computeContact( const ForceType& force, ParticleType& particles, + const ParallelType& neigh_op_tag ) +{ + auto n_local = particles.n_local; + auto f = particles.sliceForce(); + auto f_a = particles.sliceForceAtomic(); + + // Do NOT reset force - this is assumed to be done by the primary force + // kernel. + + // if ( half_neigh ) + // Forces must be atomic for half list + // compute_force_half( f_a, x, u, neigh_list, n_local, + // neigh_op_tag ); + + // Forces only atomic if using team threading + if ( std::is_same::value ) + force.computeContactFull( f_a, particles, n_local, neigh_op_tag ); + else + force.computeContactFull( f, particles, n_local, neigh_op_tag ); + Kokkos::fence(); +} + +} // namespace CabanaPD + +#endif diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index c1bdb472..9935a37e 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -228,6 +228,7 @@ class Particles _init_timer.stop(); } + // FIXME: extra option here only meant for quick route to disk impact. template void createParticles( const ExecSpace& exec_space ) { @@ -397,9 +398,9 @@ class Particles auto getForce() { return _plist_f; } auto getReferencePosition() { return _plist_x; } - void updateCurrentPosition() + void updateCurrentPosition() const { - _timer.start(); + //_timer.start(); // Not using slice function because this is called inside. auto y = Cabana::slice<0>( _aosoa_y, "current_positions" ); auto x = sliceReferencePosition(); @@ -412,7 +413,7 @@ class Particles }; Kokkos::parallel_for( "CabanaPD::CalculateCurrentPositions", policy, sum_x_u ); - _timer.stop(); + //_timer.stop(); } void resize( int new_local, int new_ghost ) diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 04329965..1750d4ff 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -669,6 +669,131 @@ class SolverFracture using base_type::print; }; +template +class SolverContact + : public SolverFracture +{ + public: + using base_type = + SolverFracture; + using exec_space = typename base_type::exec_space; + using memory_space = typename base_type::memory_space; + + using particle_type = typename base_type::particle_type; + using integrator_type = typename base_type::integrator_type; + using comm_type = typename base_type::comm_type; + 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 input_type = typename base_type::input_type; + using contact_type = Contact; + using contact_model_type = ContactModel; + + SolverContact( input_type _inputs, + std::shared_ptr _particles, + force_model_type force_model, + contact_model_type _contact_model ) + : base_type( _inputs, _particles, force_model ) + { + contact = std::make_shared( _inputs["half_neigh"], + _contact_model, *particles ); + } + + void init( const bool initial_output = true ) + { + base_type::init( false ); + + // Compute initial contact + computeContact( *contact, *particles, neigh_iter_tag{} ); + + if ( initial_output ) + particles->output( 0, 0.0, output_reference ); + } + + template + void init( BoundaryType boundary_condition, + const bool initial_output = true ) + { + base_type::init( boundary_condition, false ); + + // Compute initial contact + // FIXME: this should be before the final BC + computeContact( *contact, *particles, neigh_iter_tag{} ); + + if ( initial_output ) + particles->output( 0, 0.0, output_reference ); + } + + template + void run( BoundaryType boundary_condition ) + { + this->init_output( boundary_condition.timeInit() ); + + // Main timestep loop. + for ( int step = 1; step <= num_steps; step++ ) + { + _step_timer.start(); + + // Integrate - velocity Verlet first half. + integrator->initialHalfStep( *particles ); + + // Add non-force boundary condition. + if ( !boundary_condition.forceUpdate() ) + boundary_condition.apply( exec_space(), *particles, step * dt ); + + if constexpr ( std::is_same::value ) + comm->gatherTemperature(); + + // Update ghost particles. + comm->gatherDisplacement(); + + // Compute internal forces. + base_type::updateForce(); + + computeContact( *contact, *particles, neigh_iter_tag{} ); + + // Add force boundary condition. + if ( boundary_condition.forceUpdate() ) + boundary_condition.apply( exec_space{}, *particles, step * dt ); + + // Integrate - velocity Verlet second half. + integrator->finalHalfStep( *particles ); + + base_type::output( step ); + } + + // Final output and timings. + this->final_output(); + } + + using base_type::dt; + using base_type::num_steps; + using base_type::output_frequency; + using base_type::output_reference; + + protected: + using base_type::comm; + using base_type::force; + using base_type::inputs; + using base_type::integrator; + using base_type::neighbors; + using base_type::particles; + + using base_type::mu; + + std::shared_ptr contact; + + using base_type::_init_time; + using base_type::_init_timer; + using base_type::_step_timer; + using base_type::print; +}; + +// =============================================================== + template auto createSolverElastic( InputsType inputs, @@ -702,6 +827,17 @@ auto createSolverFracture( InputsType inputs, inputs, particles, model, prenotch ); } +template +auto createSolverContact( InputType inputs, + std::shared_ptr particles, + ForceModel model, ContactModel contact ) +{ + return std::make_shared>( + inputs, particles, model, contact ); +} + } // namespace CabanaPD #endif