From 575a524b3092c484d709cfe6d5b41283140f3c4f Mon Sep 17 00:00:00 2001 From: Austin Isner Date: Fri, 6 Dec 2024 09:16:34 -0500 Subject: [PATCH 1/3] Add HertzianContact with normal spring and damping forces; includes new test with physically realistic values for 316L --- src/CabanaPD.hpp | 2 + src/force/CabanaPD_Force_Contact.hpp | 17 ++ src/force/CabanaPD_Force_HertzianContact.hpp | 164 +++++++++++++++++++ src/force/CabanaPD_HertzianContact.hpp | 70 ++++++++ unit_test/CMakeLists.txt | 7 +- unit_test/inputs/hertzian_contact.json | 16 ++ unit_test/tstHertz.hpp | 141 ++++++++++++++++ 7 files changed, 416 insertions(+), 1 deletion(-) create mode 100644 src/force/CabanaPD_Force_HertzianContact.hpp create mode 100644 src/force/CabanaPD_HertzianContact.hpp create mode 100644 unit_test/inputs/hertzian_contact.json create mode 100644 unit_test/tstHertz.hpp diff --git a/src/CabanaPD.hpp b/src/CabanaPD.hpp index f1014885..ec810cc4 100644 --- a/src/CabanaPD.hpp +++ b/src/CabanaPD.hpp @@ -33,7 +33,9 @@ #include #include #include +#include #include #include +#include #endif diff --git a/src/force/CabanaPD_Force_Contact.hpp b/src/force/CabanaPD_Force_Contact.hpp index f3b389ef..46dd9f69 100644 --- a/src/force/CabanaPD_Force_Contact.hpp +++ b/src/force/CabanaPD_Force_Contact.hpp @@ -20,6 +20,23 @@ namespace CabanaPD { +/****************************************************************************** +Contact helper functions +******************************************************************************/ +template +KOKKOS_INLINE_FUNCTION void getRelativeNormalVelocityComponents( + const VelType& vel, const int i, const int j, const double rx, + const double ry, const double rz, const double r, double& vx, double& vy, + double& vz, double& vn ) +{ + vx = vel( i, 0 ) - vel( j, 0 ); + vy = vel( i, 1 ) - vel( j, 1 ); + vz = vel( i, 2 ) - vel( j, 2 ); + + vn = vx * rx + vy * ry + vz * rz; + vn /= r; +}; + /****************************************************************************** Normal repulsion forces ******************************************************************************/ diff --git a/src/force/CabanaPD_Force_HertzianContact.hpp b/src/force/CabanaPD_Force_HertzianContact.hpp new file mode 100644 index 00000000..30d52197 --- /dev/null +++ b/src/force/CabanaPD_Force_HertzianContact.hpp @@ -0,0 +1,164 @@ +/**************************************************************************** + * 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_HERTZIAN_H +#define CONTACT_HERTZIAN_H + +#include + +#include +#include +#include +#include +#include + +namespace CabanaPD +{ +/****************************************************************************** + Normal repulsion forces +******************************************************************************/ +template +class Force + : public Force +{ + public: + using base_type = Force; + using neighbor_list_type = typename base_type::neighbor_list_type; + + template + Force( const bool half_neigh, const ParticleType& particles, + const HertzianModel model ) + : base_type( half_neigh, model.Rc, particles.sliceCurrentPosition(), + particles.n_local, particles.ghost_mesh_lo, + particles.ghost_mesh_hi ) + , _model( model ) + { + for ( int d = 0; d < particles.dim; d++ ) + { + mesh_min[d] = particles.ghost_mesh_lo[d]; + mesh_max[d] = particles.ghost_mesh_hi[d]; + } + } + + template + void computeForceFull( ForceType& fc, const PosType& x, const PosType& u, + const ParticleType& particles, const int n_local, + ParallelType& neigh_op_tag ) + { + auto Rc = _model.Rc; + auto radius = _model.radius; + auto Es = _model.Es; + auto Rs = _model.Rs; + auto beta = _model.beta; + + const double coeff_h_n = 4.0 / 3.0 * Es * std::sqrt( Rs ); + const double coeff_h_d = -2.0 * sqrt( 5.0 / 6.0 ) * beta; + + const auto vol = particles.sliceVolume(); + const auto rho = particles.sliceDensity(); + const auto y = particles.sliceCurrentPosition(); + const auto vel = particles.sliceVelocity(); + + _neigh_timer.start(); + _neigh_list.build( y, 0, n_local, Rc, 1.0, mesh_min, mesh_max ); + _neigh_timer.stop(); + + auto contact_full = KOKKOS_LAMBDA( const int i, const int j ) + { + using Kokkos::abs; + using Kokkos::min; + using Kokkos::pow; + using Kokkos::sqrt; + + 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 "overlap" + const double delta_n = ( r - 2.0 * radius ); + + // Hertz normal force coefficient + double coeff = 0.0; + double Sn = 0.0; + if ( delta_n < 0.0 ) + { + coeff = + min( 0.0, -coeff_h_n * pow( abs( delta_n ), 3.0 / 2.0 ) ); + Sn = 2.0 * Es * sqrt( Rs * abs( delta_n ) ); + } + + coeff /= vol( i ); + + 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; + + // Hertz normal force damping component + double vx, vy, vz, vn; + getRelativeNormalVelocityComponents( vel, i, j, rx, ry, rz, r, vx, + vy, vz, vn ); + + double ms = ( rho( i ) * vol( i ) ) / 2.0; + double fnd = coeff_h_d * sqrt( Sn * ms ) * vn / vol( i ); + + fcx_i = fnd * rx / r; + fcy_i = fnd * ry / r; + fcz_i = fnd * rz / r; + + fc( i, 0 ) += fcx_i; + fc( i, 1 ) += fcy_i; + fc( i, 2 ) += fcz_i; + }; + + _timer.start(); + + // 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, _neigh_list, Cabana::FirstNeighborsTag(), + neigh_op_tag, "CabanaPD::Contact::compute_full" ); + + _timer.stop(); + } + + // FIXME: implement energy + template + double computeEnergyFull( WType&, const PosType&, const PosType&, + ParticleType&, const int, ParallelType& ) + { + return 0.0; + } + + protected: + HertzianModel _model; + using base_type::_half_neigh; + using base_type::_neigh_list; + using base_type::_timer; + Timer _neigh_timer; + + double mesh_max[3]; + double mesh_min[3]; +}; + +} // namespace CabanaPD + +#endif diff --git a/src/force/CabanaPD_HertzianContact.hpp b/src/force/CabanaPD_HertzianContact.hpp new file mode 100644 index 00000000..56435da7 --- /dev/null +++ b/src/force/CabanaPD_HertzianContact.hpp @@ -0,0 +1,70 @@ +/**************************************************************************** + * 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 CONTACTMODEL_HERTZIAN_H +#define CONTACTMODEL_HERTZIAN_H + +#include + +#include +#include +#include + +namespace CabanaPD +{ +struct HertzianModel : public ContactModel +{ + // FIXME: This is for use as the primary force model. + using base_model = PMB; + using fracture_type = Elastic; + using thermal_type = TemperatureIndependent; + + using ContactModel::Rc; // Contact horizon (should be > 2*radius) + + double nu; // Poisson's ratio + double radius; // Actual radius + double Rs; // Equivalent radius + double Es; // Equivalent Young's modulus + double e; // Coefficient of restitution + double beta; // Damping coefficient + + HertzianModel( const double _Rc, const double _radius, const double _nu, + const double _E, const double _e ) + : ContactModel( 1.0, _Rc ) + { + set_param( _radius, _nu, _E, _e ); + } + + void set_param( const double _radius, const double _nu, const double _E, + const double _e ) + { + using Kokkos::log; + using Kokkos::pow; + using Kokkos::sqrt; + + nu = _nu; + radius = _radius; + Rs = 0.5 * radius; + Es = _E / ( 2.0 * pow( 1.0 - nu, 2.0 ) ); + e = _e; + double ln_e = log( e ); + beta = -ln_e / sqrt( pow( ln_e, 2.0 ) + pow( pi, 2.0 ) ); + } +}; + +template <> +struct is_contact : public std::true_type +{ +}; + +} // namespace CabanaPD + +#endif diff --git a/unit_test/CMakeLists.txt b/unit_test/CMakeLists.txt index b97153e3..d3003f14 100644 --- a/unit_test/CMakeLists.txt +++ b/unit_test/CMakeLists.txt @@ -1,3 +1,8 @@ +configure_file( + ${CMAKE_CURRENT_SOURCE_DIR}/inputs/hertzian_contact.json + ${CMAKE_CURRENT_BINARY_DIR}/hertzian_contact.json + COPYONLY +) ##--------------------------------------------------------------------------## ## On-node tests ##--------------------------------------------------------------------------## @@ -48,6 +53,6 @@ macro(CabanaPD_add_tests) endforeach() endmacro() -CabanaPD_add_tests(NAMES Force Integrator) +CabanaPD_add_tests(NAMES Force Integrator Hertz) CabanaPD_add_tests(MPI NAMES Comm) diff --git a/unit_test/inputs/hertzian_contact.json b/unit_test/inputs/hertzian_contact.json new file mode 100644 index 00000000..5180c684 --- /dev/null +++ b/unit_test/inputs/hertzian_contact.json @@ -0,0 +1,16 @@ +{ + "num_cells" : {"value": [1, 1, 1]}, + "system_size" : {"value": [1.0e-3, 1.0e-3, 1.0e-3], "unit": "m"}, + "density" : {"value": 7.95e3, "unit": "kg/m^3"}, + "volume" : {"value": 5.236e-13, "unit": "m^3"}, + "elastic_modulus" : {"value": 195.6e9, "unit": "Pa"}, + "poisson_ratio" : {"value": 0.25, "unit": ""}, + "restitution" : {"value": 0.5}, + "horizon" : {"value": 1e-4 }, + "radius" : {"value": 5e-5 }, + "final_time" : {"value": 1e-5, "unit": "s"}, + "timestep" : {"value": 1e-8, "unit": "s"}, + "timestep_safety_factor" : {"value": 0.8}, + "output_frequency" : {"value": 10000}, + "output_reference" : {"value": false} +} diff --git a/unit_test/tstHertz.hpp b/unit_test/tstHertz.hpp new file mode 100644 index 00000000..f9e08bb3 --- /dev/null +++ b/unit_test/tstHertz.hpp @@ -0,0 +1,141 @@ +/**************************************************************************** + * 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 * + ****************************************************************************/ + +#include + +#include +#include + +namespace Test +{ +template +double calculateKE( const VelType& v, const DensityType& rho, + const VolumeType& vol ) +{ + using Kokkos::hypot; + using Kokkos::pow; + + double tke = 0.0; + Kokkos::parallel_reduce( + "total_ke", v.size(), + KOKKOS_LAMBDA( const int i, double& sum ) { + sum += 0.5 * rho( i ) * vol( i ) * + pow( hypot( v( i, 0 ), v( i, 1 ), v( i, 2 ) ), 2.0 ); + }, + tke ); + return tke; +} + +void testHertzianContact( const std::string filename ) +{ + // ==================================================== + // Use test Kokkos spaces + // ==================================================== + using exec_space = TEST_EXECSPACE; + using memory_space = TEST_MEMSPACE; + + // ==================================================== + // Read inputs + // ==================================================== + CabanaPD::Inputs inputs( filename ); + + // ==================================================== + // Material parameters + // ==================================================== + double rho0 = inputs["density"]; + double vol = inputs["volume"]; + double delta = inputs["horizon"]; + delta += 1e-10; + double radius = inputs["radius"]; + double nu = inputs["poisson_ratio"]; + double E = inputs["elastic_modulus"]; + double e = inputs["restitution"]; + + // ==================================================== + // Discretization + // ==================================================== + std::array low_corner = inputs["low_corner"]; + std::array high_corner = inputs["high_corner"]; + std::array num_cells = inputs["num_cells"]; + + // ==================================================== + // Custom particle creation + // ==================================================== + const int num_particles = 2; + // Purposely using zero-init here. + Kokkos::View position( "custom_position", 2 ); + Kokkos::View volume( "custom_volume", 2 ); + + Kokkos::parallel_for( + "create_particles", Kokkos::RangePolicy( 0, num_particles ), + KOKKOS_LAMBDA( const int p ) { + if ( p == 0 ) + position( p, 0 ) = 0.51e-4; + else + position( p, 0 ) = -0.51e-4; + volume( p ) = vol; + } ); + + // ==================================================== + // Force model + // ==================================================== + using model_type = CabanaPD::HertzianModel; + model_type contact_model( delta, radius, nu, E, e ); + + // ==================================================== + // Particle generation + // ==================================================== + int halo_width = 1; + auto particles = CabanaPD::createParticles( + exec_space{}, position, volume, low_corner, high_corner, num_cells, + halo_width ); + + // ==================================================== + // Custom particle initialization + // ==================================================== + auto rho = particles->sliceDensity(); + auto v = particles->sliceVelocity(); + auto vo = particles->sliceVolume(); + + auto init_functor = KOKKOS_LAMBDA( const int p ) + { + // Density + rho( p ) = rho0; + if ( p == 0 ) + v( p, 0 ) = -1.0; + else + v( p, 0 ) = 1.0; + }; + particles->updateParticles( exec_space{}, init_functor ); + + // Get initial total KE + double ke_i = calculateKE( v, rho, vo ); + + // ==================================================== + // Simulation run + // ==================================================== + auto cabana_pd = CabanaPD::createSolverElastic( + inputs, particles, contact_model ); + cabana_pd->init(); + cabana_pd->run(); + + // Get final total KE + double ke_f = calculateKE( v, rho, vo ); + EXPECT_NEAR( std::sqrt( ke_f / ke_i ), e, 1e-3 ); +} + +TEST( TEST_CATEGORY, test_hertzian_contact ) +{ + std::string input = "hertzian_contact.json"; + testHertzianContact( input ); +} + +} // end namespace Test From 5bf465190bb8b47c0b3fc4c928bec4ea4effc157 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 9 Dec 2024 18:01:26 -0500 Subject: [PATCH 2/3] Add custom particle position/volume creation --- src/CabanaPD_Particles.hpp | 113 ++++++++++++++++++++++++++++++++++--- 1 file changed, 106 insertions(+), 7 deletions(-) diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 84b87323..a85541c3 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -188,6 +188,23 @@ class Particles createParticles( exec_space, user_create ); } + // Constructor with existing particle data. + template + Particles( const ExecSpace& exec_space, const PositionType& x, + const VolumeType& vol, std::array low_corner, + std::array high_corner, + const std::array num_cells, const int max_halo_width ) + : halo_width( max_halo_width ) + , _plist_x( "positions" ) + , _plist_f( "forces" ) + { + createDomain( low_corner, high_corner, num_cells ); + + _init_timer.start(); + createParticles( exec_space, x, vol ); + _init_timer.stop(); + } + void createDomain( std::array low_corner, std::array high_corner, const std::array num_cells ) @@ -306,6 +323,44 @@ class Particles _init_timer.stop(); } + // Store custom created particle positions and volumes. + template + void createParticles( const ExecSpace, const PositionType& x, + const VolumeType& vol ) + { + resize( vol.size(), 0 ); + auto p_x = sliceReferencePosition(); + auto p_vol = sliceVolume(); + auto v = sliceVelocity(); + auto f = sliceForce(); + auto type = sliceType(); + auto rho = sliceDensity(); + auto u = sliceDisplacement(); + auto nofail = sliceNoFail(); + + static_assert( + Cabana::is_accessible_from< + memory_space, typename PositionType::execution_space>::value ); + + Kokkos::parallel_for( + "copy_to_particles", Kokkos::RangePolicy( 0, n_local ), + KOKKOS_LAMBDA( const int pid ) { + // Set the particle position and volume. + // Set everything else to zero. + p_vol( pid ) = vol( pid ); + for ( int d = 0; d < 3; d++ ) + { + p_x( pid, d ) = x( pid, d ); + u( pid, d ) = 0.0; + v( pid, d ) = 0.0; + f( pid, d ) = 0.0; + } + type( pid ) = 0; + nofail( pid ) = 0; + rho( pid ) = 1.0; + } ); + } + template void updateParticles( const ExecSpace, const FunctorType init_functor ) { @@ -798,13 +853,10 @@ class Particles _aosoa_output = aosoa_output_type( "Particle Output Fields", 0 ); } - // Constructor which initializes particles on regular grid. - template - Particles( const ExecSpace& exec_space, std::array low_corner, - std::array high_corner, - const std::array num_cells, const int max_halo_width ) - : base_type( exec_space, low_corner, high_corner, num_cells, - max_halo_width ) + // Base constructor. + template + Particles( Args&&... args ) + : base_type( std::forward( args )... ) { _aosoa_output = aosoa_output_type( "Particle Output Fields", n_local ); init_output(); @@ -945,6 +997,53 @@ auto createParticles( exec_space, low_corner, high_corner, num_cells, max_halo_width, EnergyOutput{} ); } + +template +auto createParticles( const ExecSpace& exec_space, const PositionType& x, + const VolumeType& vol, std::array low_corner, + std::array high_corner, + const std::array num_cells, + const int max_halo_width, OutputType ) +{ + return std::make_shared< + CabanaPD::Particles>( + exec_space, x, vol, low_corner, high_corner, num_cells, + max_halo_width ); +} + +template +auto createParticles( + const ExecSpace& exec_space, const PositionType& x, const VolumeType& vol, + std::array low_corner, std::array high_corner, + const std::array num_cells, const int max_halo_width, OutputType, + typename std::enable_if<( is_temperature_dependent::value ), + int>::type* = 0 ) +{ + return std::make_shared>( + exec_space, x, vol, low_corner, high_corner, num_cells, + max_halo_width ); +} + +template +auto createParticles( const ExecSpace& exec_space, const PositionType& x, + const VolumeType& vol, std::array low_corner, + std::array high_corner, + const std::array num_cells, + const int max_halo_width ) +{ + return createParticles( exec_space, x, vol, low_corner, + high_corner, num_cells, + max_halo_width, EnergyOutput{} ); +} + } // namespace CabanaPD #endif From 35c1d128d4ec2122a53183b72a0419da92b7c030 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 18 Dec 2024 18:25:56 -0500 Subject: [PATCH 3/3] Bump min Kokkos to 4.1 for math functions; bump min Cabana to 0.7 for simplicity --- .github/workflows/CI.yml | 22 ++++++++++------------ CMakeLists.txt | 2 +- 2 files changed, 11 insertions(+), 13 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 01f68736..c792677e 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -20,26 +20,26 @@ jobs: cxx: ['g++', 'clang++'] backend: ['SERIAL', 'OPENMP'] cmake_build_type: ['Debug', 'Release'] - kokkos_ver: ['3.7.02'] + kokkos_ver: ['4.1.00'] output: ['HDF5'] include: - distro: 'ubuntu:latest' cxx: 'g++' backend: 'SERIAL' cmake_build_type: 'Debug' - kokkos_ver: '3.7.02' + kokkos_ver: '4.1.00' output: 'SILO' - distro: 'ubuntu:latest' cxx: 'g++' backend: 'SERIAL' cmake_build_type: 'Debug' - kokkos_ver: '3.7.02' + kokkos_ver: '4.1.00' output: 'NONE' - distro: 'ubuntu:latest' cxx: 'g++' backend: 'SERIAL' cmake_build_type: 'Debug' - kokkos_ver: '3.7.02' + kokkos_ver: '4.1.00' output: 'BOTH' runs-on: ubuntu-20.04 container: ghcr.io/ecp-copa/ci-containers/${{ matrix.distro }} @@ -65,8 +65,7 @@ jobs: uses: actions/checkout@v3 with: repository: ECP-CoPA/Cabana - # This version is post-release 0.6 - ref: f99c7db9d54c57373ada6b16132c20d89d1ebb8e + ref: 0.7.0 path: cabana - name: Build Cabana working-directory: cabana @@ -118,7 +117,8 @@ jobs: matrix: cxx: ['hipcc'] cmake_build_type: ['Release'] - kokkos_ver: ['3.7.02'] + # Using >4.1 because of kokkos build error without available device + kokkos_ver: ['4.2.01'] runs-on: ubuntu-20.04 container: ghcr.io/ecp-copa/ci-containers/rocm:latest steps: @@ -158,8 +158,7 @@ jobs: uses: actions/checkout@v3 with: repository: ECP-CoPA/Cabana - # This version is post-release 0.6 - ref: f99c7db9d54c57373ada6b16132c20d89d1ebb8e + ref: 0.7.0 path: cabana - name: Build Cabana working-directory: cabana @@ -200,7 +199,7 @@ jobs: strategy: matrix: cmake_build_type: ['Release'] - kokkos_ver: ['3.7.02'] + kokkos_ver: ['4.1.00'] runs-on: ubuntu-20.04 container: ghcr.io/ecp-copa/ci-containers/cuda:12.2.0 steps: @@ -240,8 +239,7 @@ jobs: uses: actions/checkout@v3 with: repository: ECP-CoPA/Cabana - # This version is post-release 0.6 - ref: f99c7db9d54c57373ada6b16132c20d89d1ebb8e + ref: 0.7.0 path: cabana - name: Build Cabana working-directory: cabana diff --git a/CMakeLists.txt b/CMakeLists.txt index ec5c4337..a0e62fbe 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -17,7 +17,7 @@ set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake) ##---------------------------------------------------------------------------## # Set up main options (inherit from Kokkos and Cabana CMake) ##---------------------------------------------------------------------------## -find_package(Cabana REQUIRED) +find_package(Cabana REQUIRED 0.7.0) macro(CabanaPD_check_optional) cmake_parse_arguments(CABANA "" "OPTION" "" ${ARGN})