Initial contact forces

Co-authored-by: Pablo Seleson <[email protected]>

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})

examples/mechanics/disk_impact.cpp

@@ -0,0 +1,101 @@
+#include <cmath>
+#include <fstream>
+#include <iostream>
+
+#include "mpi.h"
+
+#include <Kokkos_Core.hpp>
+
+#include <CabanaPD.hpp>
+
+// 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<CabanaPD::PMB, CabanaPD::Fracture>;
+    model_type force_model( delta, K, G0 );
+
+    // Create particles from mesh.
+    auto particles = CabanaPD::createParticles<memory_space, model_type>(
+        exec_space(), inputs );
+
+    std::array<double, 3> low_corner = inputs["low_corner"];
+    std::array<double, 3> 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<memory_space>(
+        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();
+}

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"}
+}

src/CabanaPD.hpp

@@ -16,6 +16,7 @@
 #include <CabanaPD_Boundary.hpp>
 #include <CabanaPD_Comm.hpp>
 #include <CabanaPD_Constants.hpp>
+#include <CabanaPD_Contact.hpp>
 #include <CabanaPD_DisplacementProfile.hpp>
 #include <CabanaPD_Fields.hpp>
 #include <CabanaPD_Force.hpp>

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 <cmath>
+
+#include <CabanaPD_Force.hpp>
+#include <CabanaPD_Input.hpp>
+#include <CabanaPD_Output.hpp>
+
+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 MemorySpace, class ContactType>
+class Contact;
+
+/******************************************************************************
+  Normal repulsion computation
+******************************************************************************/
+template <class MemorySpace>
+class Contact<MemorySpace, NormalRepulsionModel>
+{
+  public:
+    using neighbor_type =
+        Cabana::VerletList<MemorySpace, Cabana::FullNeighborTag,
+                           Cabana::VerletLayout2D, Cabana::TeamOpTag>;
+
+    template <class ParticleType>
+    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<neighbor_type>(
+            y, 0, particles.n_local, model.Rc, 1.0, mesh_min, mesh_max );
+    }
+
+    template <class ForceType, class ParticleType, class ParallelType>
+    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<exec_space> 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<neighbor_type> _neighbors;
+
+    double mesh_max[3];
+    double mesh_min[3];
+};
+
+template <class ForceType, class ParticleType, class ParallelType>
+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<decltype( neigh_op_tag ), Cabana::TeamOpTag>::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

src/CabanaPD_Particles.hpp

@@ -228,6 +228,7 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, Dimension>
         _init_timer.stop();
     }
 
+    // FIXME: extra option here only meant for quick route to disk impact.
     template <class ExecSpace>
     void createParticles( const ExecSpace& exec_space )
     {
@@ -397,9 +398,9 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, Dimension>
     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<MemorySpace, PMB, TemperatureIndependent, Dimension>
         };
         Kokkos::parallel_for( "CabanaPD::CalculateCurrentPositions", policy,
                               sum_x_u );
-        _timer.stop();
+        //_timer.stop();
     }
 
     void resize( int new_local, int new_ghost )