Merge pull request #158 from abisner/dem-model

DEM Force Model

.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

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

src/CabanaPD.hpp

@@ -33,7 +33,9 @@
 #include <force/CabanaPD_ForceModels_LPS.hpp>
 #include <force/CabanaPD_ForceModels_PMB.hpp>
 #include <force/CabanaPD_Force_Contact.hpp>
+#include <force/CabanaPD_Force_HertzianContact.hpp>
 #include <force/CabanaPD_Force_LPS.hpp>
 #include <force/CabanaPD_Force_PMB.hpp>
+#include <force/CabanaPD_HertzianContact.hpp>
 
 #endif

src/CabanaPD_Particles.hpp

@@ -188,6 +188,23 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, BaseOutput, Dimension>
         createParticles( exec_space, user_create );
     }
 
+    // Constructor with existing particle data.
+    template <class ExecSpace, class PositionType, class VolumeType>
+    Particles( const ExecSpace& exec_space, const PositionType& x,
+               const VolumeType& vol, std::array<double, dim> low_corner,
+               std::array<double, dim> high_corner,
+               const std::array<int, dim> 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<double, dim> low_corner,
                        std::array<double, dim> high_corner,
                        const std::array<int, dim> num_cells )
@@ -306,6 +323,44 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, BaseOutput, Dimension>
         _init_timer.stop();
     }
 
+    // Store custom created particle positions and volumes.
+    template <class ExecSpace, class PositionType, class VolumeType>
+    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<ExecSpace>( 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 <class ExecSpace, class FunctorType>
     void updateParticles( const ExecSpace, const FunctorType init_functor )
     {
@@ -798,13 +853,10 @@ class Particles<MemorySpace, ModelType, ThermalType, EnergyOutput, Dimension>
         _aosoa_output = aosoa_output_type( "Particle Output Fields", 0 );
     }
 
-    // Constructor which initializes particles on regular grid.
-    template <class ExecSpace>
-    Particles( const ExecSpace& exec_space, std::array<double, dim> low_corner,
-               std::array<double, dim> high_corner,
-               const std::array<int, dim> num_cells, const int max_halo_width )
-        : base_type( exec_space, low_corner, high_corner, num_cells,
-                     max_halo_width )
+    // Base constructor.
+    template <typename... Args>
+    Particles( Args&&... args )
+        : base_type( std::forward<Args>( 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 <typename MemorySpace, typename ModelType, typename ExecSpace,
+          typename PositionType, typename VolumeType, std::size_t Dim,
+          typename OutputType>
+auto createParticles( const ExecSpace& exec_space, const PositionType& x,
+                      const VolumeType& vol, std::array<double, Dim> low_corner,
+                      std::array<double, Dim> high_corner,
+                      const std::array<int, Dim> num_cells,
+                      const int max_halo_width, OutputType )
+{
+    return std::make_shared<
+        CabanaPD::Particles<MemorySpace, typename ModelType::base_model,
+                            typename ModelType::thermal_type, OutputType>>(
+        exec_space, x, vol, low_corner, high_corner, num_cells,
+        max_halo_width );
+}
+
+template <typename MemorySpace, typename ModelType, typename ThermalType,
+          typename ExecSpace, typename PositionType, typename VolumeType,
+          std::size_t Dim, typename OutputType>
+auto createParticles(
+    const ExecSpace& exec_space, const PositionType& x, const VolumeType& vol,
+    std::array<double, Dim> low_corner, std::array<double, Dim> high_corner,
+    const std::array<int, Dim> num_cells, const int max_halo_width, OutputType,
+    typename std::enable_if<( is_temperature_dependent<ThermalType>::value ),
+                            int>::type* = 0 )
+{
+    return std::make_shared<CabanaPD::Particles<
+        MemorySpace, ModelType, typename ThermalType::base_type, OutputType>>(
+        exec_space, x, vol, low_corner, high_corner, num_cells,
+        max_halo_width );
+}
+
+template <typename MemorySpace, typename ModelType, typename ExecSpace,
+          typename PositionType, typename VolumeType, std::size_t Dim>
+auto createParticles( const ExecSpace& exec_space, const PositionType& x,
+                      const VolumeType& vol, std::array<double, Dim> low_corner,
+                      std::array<double, Dim> high_corner,
+                      const std::array<int, Dim> num_cells,
+                      const int max_halo_width )
+{
+    return createParticles<MemorySpace, ModelType, ExecSpace, PositionType,
+                           VolumeType, Dim>( exec_space, x, vol, low_corner,
+                                             high_corner, num_cells,
+                                             max_halo_width, EnergyOutput{} );
+}
+
 } // namespace CabanaPD
 
 #endif

src/force/CabanaPD_Force_Contact.hpp

@@ -20,6 +20,23 @@
 
 namespace CabanaPD
 {
+/******************************************************************************
+Contact helper functions
+******************************************************************************/
+template <class VelType>
+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
 ******************************************************************************/