Impactor example

examples/CMakeLists.txt

@@ -7,4 +7,7 @@ target_link_libraries(KalthoffWinkler LINK_PUBLIC CabanaPD)
 add_executable(CrackBranching crack_branching.cpp)
 target_link_libraries(CrackBranching LINK_PUBLIC CabanaPD)
 
-install(TARGETS ElasticWave KalthoffWinkler CrackBranching DESTINATION ${CMAKE_INSTALL_BINDIR})
+add_executable(DiskImpact disk_impact.cpp)
+target_link_libraries(DiskImpact LINK_PUBLIC CabanaPD)
+
+install(TARGETS ElasticWave KalthoffWinkler CrackBranching DiskImpact DESTINATION ${CMAKE_INSTALL_BINDIR})

examples/disk_impact.cpp

@@ -0,0 +1,91 @@
+#include <cmath>
+#include <fstream>
+#include <iostream>
+
+#include "mpi.h"
+
+#include <Kokkos_Core.hpp>
+
+#include <CabanaPD.hpp>
+
+int main( int argc, char* argv[] )
+{
+    MPI_Init( &argc, &argv );
+
+    Kokkos::initialize( argc, argv );
+
+    // FIXME: change backend at compile time for now.
+    using exec_space = Kokkos::DefaultExecutionSpace;
+    using memory_space = typename exec_space::memory_space;
+
+    // Plate dimension)
+    double radius = 37.0e-3;   // [m]
+    double thickness = 2.5e-3; // [m]
+
+    // Domain
+    std::array<int, 3> num_cell = { 149, 149, 5 };
+    std::array<double, 3> low_corner = { -radius, -radius, -thickness };
+    std::array<double, 3> high_corner = { radius, radius, 0.0 };
+    double t_final = 2.0e-4;
+    double dt = 1.0e-7;
+    int output_frequency = 50;
+    bool output_reference = false;
+
+    // Material constants
+    // double nu = 1.0 / 4.0; // unitless
+    double K = 14.9e+9;   // [Pa]
+    double rho0 = 2200.0; // [kg/m^3]
+    double G0 = 10.0575;  // [J/m^2]
+    double impact_v = -100.0;
+    // double ks = 1.0e17;
+    double impact_r = 5e-3; // [m]
+
+    // FIXME: set halo width based on delta
+    double delta = 0.0015;
+    int halo_width = 2;
+
+    // Choose force model type.
+    using model_type = CabanaPD::ForceModel<CabanaPD::PMB, CabanaPD::Fracture>;
+    model_type force_model( delta, K, G0 );
+    CabanaPD::Inputs<3> inputs( num_cell, low_corner, high_corner, t_final, dt,
+                                output_frequency, output_reference );
+    inputs.read_args( argc, argv );
+
+    // Create particles from mesh.
+    // Does not set displacements, velocities, etc.
+    auto particles = std::make_shared<
+        CabanaPD::Particles<memory_space, typename model_type::base_model>>(
+        exec_space(), inputs.low_corner, inputs.high_corner, inputs.num_cells,
+        halo_width, true );
+    double dx = particles->dx[0];
+    double r_c = dx * 2.0;
+
+    // Define particle initialization.
+    auto x = particles->sliceReferencePosition();
+    auto v = particles->sliceVelocity();
+    auto f = particles->sliceForce();
+    auto rho = particles->sliceDensity();
+
+    CabanaPD::Prenotch<0> prenotch;
+    CabanaPD::BoundaryCondition<CabanaPD::AllTag, CabanaPD::ImpactBCTag> bc(
+        impact_r, impact_v, dt );
+
+    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 );
+
+    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, bc, prenotch, contact_model );
+    cabana_pd->init_force();
+    cabana_pd->run();
+
+    Kokkos::finalize();
+    MPI_Finalize();
+}

src/CabanaPD.hpp

@@ -14,6 +14,7 @@
 
 #include <CabanaPD_Boundary.hpp>
 #include <CabanaPD_Comm.hpp>
+#include <CabanaPD_Contact.hpp>
 #include <CabanaPD_Force.hpp>
 #include <CabanaPD_Input.hpp>
 #include <CabanaPD_Integrate.hpp>

src/CabanaPD_Boundary.hpp

@@ -262,6 +262,59 @@ auto createBoundaryCondition( BCTag, ExecSpace exec_space, Particles particles,
     return BoundaryCondition<bc_index_type, BCTag>( value, bc_indices );
 }
 
+struct ImpactBCTag
+{
+};
+struct AllTag
+{
+};
+
+template <class AllTag, class BCTag>
+struct BoundaryCondition;
+
+template <>
+struct BoundaryCondition<AllTag, ImpactBCTag>
+{
+    double _R;
+    double _v;
+    double _dt;
+    double _x;
+    double _y;
+    double _z;
+
+    BoundaryCondition( const double R, const double v, const double dt )
+        : _R( R )
+        , _v( v )
+        , _dt( dt )
+        , _x( 0.0 )
+        , _y( 0.0 )
+        , _z( R )
+    {
+    }
+
+    template <class ExecSpace, class ParticleType>
+    void apply( ExecSpace, ParticleType& particles )
+    {
+        _z += _dt * _v;
+        auto f = particles.sliceForce();
+        auto x = particles.sliceReferencePosition();
+        Kokkos::RangePolicy<ExecSpace> policy( 0, x.size() );
+        Kokkos::parallel_for(
+            "CabanaPD::BC::apply", policy, KOKKOS_LAMBDA( const int p ) {
+                double r = sqrt( ( x( p, 0 ) - _x ) * ( x( p, 0 ) - _x ) +
+                                 ( x( p, 1 ) - _y ) * ( x( p, 1 ) - _y ) +
+                                 ( x( p, 2 ) - _z ) * ( x( p, 2 ) - _z ) );
+                if ( r < _R )
+                {
+                    double fmag = -1.0e17 * ( r - _R ) * ( r - _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;
+                }
+            } );
+    }
+};
+
 } // namespace CabanaPD
 
 #endif

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