Version of contact that compiles

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,92 @@
+#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 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.
+    using device_type = Kokkos::Device<exec_space, memory_space>;
+    auto particles = std::make_shared<
+        CabanaPD::Particles<device_type, typename model_type::base_model>>(
+        exec_space(), inputs.low_corner, inputs.high_corner, inputs.num_cells,
+        halo_width, true );
+    double dx = particles->dx;
+    double r_c = dx * 2.0;
+
+    // Define particle initialization.
+    auto x = particles->sliceRefPosition();
+    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<device_type>(
+        inputs, particles, force_model, bc, prenotch, contact_model );
+    cabana_pd->init_force();
+    cabana_pd->run();
+
+    Kokkos::finalize();
+    MPI_Finalize();
+}

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.sliceRefPosition();
+        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

@@ -14,9 +14,9 @@
 
 #include <cmath>
 
+#include <CabanaPD_Force.hpp>
 #include <CabanaPD_Input.hpp>
 #include <CabanaPD_Output.hpp>
-#include <CabanaPD_Force.hpp>
 
 namespace CabanaPD
 {
@@ -29,9 +29,11 @@ struct ContactModel
     double Rc;
 
     ContactModel(){};
+    // PD horizon
+    // Contact radius
     ContactModel( const double _delta, const double _Rc )
-        : delta( _delta ){},  // PD horizon
-         Rc( _Rc ){};         // Contact radius
+        : delta( _delta )
+        , Rc( _Rc ){};
 };
 
 /* Normal repulsion */
@@ -45,8 +47,9 @@ struct NormalRepulsionModel : public ContactModel
     double K;
 
     NormalRepulsionModel(){};
-    NormalRepulsionModel( const double delta, const double Rc, const double K )
-        : ContactModel( delta ), ContactModel( Rc )
+    NormalRepulsionModel( const double delta, const double Rc, const double _K )
+        : ContactModel( delta, Rc )
+        , K( _K )
     {
         set_param( delta, Rc, K );
     }
@@ -61,54 +64,52 @@ struct NormalRepulsionModel : public ContactModel
     }
 };
 
-template <class ExecutionSpace, class ContactType>
+template <class MemorySpace, class ContactType>
 class Contact;
 
 /******************************************************************************
   Normal repulsion computation
 ******************************************************************************/
-template <class ExecutionSpace>
-class Contact<ExecutionSpace, NormalRepulsionModel>
+template <class MemorySpace>
+class Contact<MemorySpace, NormalRepulsionModel>
 {
   public:
-    using exec_space = ExecutionSpace;
     using neighbor_type =
-        Cabana::VerletList<memory_space, Cabana::FullNeighborTag,
+        Cabana::VerletList<MemorySpace, Cabana::FullNeighborTag,
                            Cabana::VerletLayout2D, Cabana::TeamOpTag>;
 
-    Contact( const bool half_neigh, const NormalRepulsionModel model )
+    template <class ParticleType>
+    Contact( const bool half_neigh, const NormalRepulsionModel model,
+             ParticleType& particles )
         : _half_neigh( half_neigh )
         , _model( model )
     {
         // Create contact neighbor list
-        double mesh_min[3] = { particles->ghost_mesh_lo[0],
-                               particles->ghost_mesh_lo[1],
-                               particles->ghost_mesh_lo[2] };
-        double mesh_max[3] = { particles->ghost_mesh_hi[0],
-                               particles->ghost_mesh_hi[1],
-                               particles->ghost_mesh_hi[2] };
-        auto x = particles->slice_x();
-        auto u = particles->slice_u();
-        auto y = particles->slice_y();
-        _contact_neighbors = std::make_shared<neighbor_type>( y, 0, particles->n_local,
-                                                     contact_model.Rc, 1.0,
-                                                     mesh_min, mesh_max );
+        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 ContactType, class PosType, class ParticleType,
-              class ParallelType>
-    void compute_contact_full( ContactType& fc, const PosType& x, const PosType& u,
-                             const ParticleType& particles,
+    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.slice_vol();
+        const auto vol = particles.sliceVolume();
+        const auto x = particles.sliceRefPosition();
+        const auto u = particles.sliceDisplacement();
+        const auto y = particles.sliceCurrentPosition();
 
-        _contact_neighbors.build( y, 0, particles->n_local,contact_model.Rc, 1.0,
-                                  mesh_min, mesh_max );
+        _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 )
         {
@@ -121,11 +122,11 @@ class Contact<ExecutionSpace, NormalRepulsionModel>
             getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz );
 
             // Contact "stretch"
-            const double sc = (r - Rc)/delta;
+            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; 
+            fcx_i = coeff * rx / r;
             fcy_i = coeff * ry / r;
             fcz_i = coeff * rz / r;
 
@@ -134,36 +135,33 @@ class Contact<ExecutionSpace, NormalRepulsionModel>
             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, _contact_neighbors, Cabana::FirstNeighborsTag(),
+            policy, contact_full, *_neighbors, Cabana::FirstNeighborsTag(),
             neigh_op_tag, "CabanaPD::Contact::compute_full" );
     }
-    
+
   protected:
     bool _half_neigh;
     NormalRepulsionModel _model;
-    neighbor_type _contact_neighbors;
+    std::shared_ptr<neighbor_type> _neighbors;
 
+    double mesh_max[3];
+    double mesh_min[3];
 };
 
-// NEED Half force computation
-
-/*
-template <class ForceType, class ParticleType, class NeighListType,
-          class ParallelType>
-void compute_force( const ForceType& force, ParticleType& particles,
-                    const NeighListType& neigh_list,
-                    const ParallelType& neigh_op_tag )
+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 x = particles.slice_x();
-    auto u = particles.slice_u();
-    auto f = particles.slice_f();
-    auto f_a = particles.slice_f_a();
+    auto f = particles.sliceForce();
+    auto f_a = particles.sliceForceAtomic();
 
-    // Reset force.
-    Cabana::deep_copy( f, 0.0 );
+    // 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
@@ -172,15 +170,12 @@ void compute_force( const ForceType& force, ParticleType& particles,
 
     // Forces only atomic if using team threading
     if ( std::is_same<decltype( neigh_op_tag ), Cabana::TeamOpTag>::value )
-        force.compute_force_full( f_a, x, u, particles, neigh_list, n_local,
-                                  neigh_op_tag );
+        force.computeContactFull( f_a, particles, n_local, neigh_op_tag );
     else
-        force.compute_force_full( f, x, u, particles, neigh_list, n_local,
-                                  neigh_op_tag );
+        force.computeContactFull( f, particles, n_local, neigh_op_tag );
     Kokkos::fence();
 }
-*/
 
 } // namespace CabanaPD
 
-#endif
+#endif