Powder fill

examples/CMakeLists.txt

@@ -1,2 +1,3 @@
 add_subdirectory(mechanics)
 add_subdirectory(thermomechanics)
+add_subdirectory(dem)

examples/dem/CMakeLists.txt

@@ -0,0 +1,3 @@
+add_executable(PowderFill powder_fill.cpp)
+target_link_libraries(PowderFill LINK_PUBLIC CabanaPD)
+install(TARGETS PowderFill DESTINATION ${CMAKE_INSTALL_BINDIR})

examples/dem/inputs/powder_fill.json

@@ -0,0 +1,19 @@
+{
+    "num_cells"                : {"value": [201, 201, 21]},
+    "system_size"              : {"value": [0.23749, 0.23749, 0.024906], "unit": "m"},
+    "density"                : {"value": 7.95e3,  "unit": "kg/m^3"},
+    "volume"                 : {"value": 5.236e-10,  "unit": "m^3"},
+    "elastic_modulus"        : {"value": 195.6e6, "unit": "Pa"},
+    "poisson_ratio"          : {"value": 0.25, "unit": ""},
+    "restitution"            : {"value": 0.5},  
+    "radius"                 : {"value": 5e-4 },
+    "horizon"                  : {"value": 0.0035, "unit": "m"},  
+    "outer_cylinder_diameter"  : {"value": 0.23749, "unit": "m"},
+    "inner_cylinder_diameter"  : {"value": 0.13589, "unit": "m"}, 
+    "wall_thickness"           : {"value": 0.0045, "unit": "m"},
+    "final_time"               : {"value": 4e0, "unit": "s"},
+    "timestep"                 : {"value": 2e-6, "unit": "s"},
+    "timestep_safety_factor"   : {"value": 0.9},
+    "output_frequency"         : {"value": 1000},
+    "output_reference"         : {"value": false}
+}

examples/dem/powder_fill.cpp

@@ -0,0 +1,179 @@
+/****************************************************************************
+ * 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 <fstream>
+#include <iostream>
+
+#include "mpi.h"
+
+#include <Kokkos_Core.hpp>
+
+#include <CabanaPD.hpp>
+
+// Simulate powder settling.
+void powderSettlingExample( const std::string filename )
+{
+    // ====================================================
+    //             Use default Kokkos spaces
+    // ====================================================
+    using exec_space = Kokkos::DefaultExecutionSpace;
+    using memory_space = typename exec_space::memory_space;
+
+    // ====================================================
+    //                   Read inputs
+    // ====================================================
+    CabanaPD::Inputs inputs( filename );
+
+    // ====================================================
+    //                Material parameters
+    // ====================================================
+    double rho0 = inputs["density"];
+    double vol0 = 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
+    // ====================================================
+    // FIXME: set halo width based on delta
+    std::array<double, 3> low_corner = inputs["low_corner"];
+    std::array<double, 3> high_corner = inputs["high_corner"];
+    std::array<int, 3> num_cells = inputs["num_cells"];
+    int m = std::floor( delta /
+                        ( ( high_corner[0] - low_corner[0] ) / num_cells[0] ) );
+    int halo_width = m + 1; // Just to be safe.
+
+    // ====================================================
+    //                    Force model
+    // ====================================================
+    using model_type = CabanaPD::HertzianModel;
+    model_type contact_model( delta / m * 0.9, radius, nu, E, e );
+
+    // ====================================================
+    //            Custom particle initialization
+    // ====================================================
+    double d_out = inputs["outer_cylinder_diameter"];
+    double d_in = inputs["inner_cylinder_diameter"];
+    double Rout = 0.5 * d_out;
+    double Rin = 0.5 * d_in;
+    double Wall_th = inputs["wall_thickness"];
+
+    // Create container.
+    auto create_container = KOKKOS_LAMBDA( const int, const double x[3] )
+    {
+        double rsq = x[0] * x[0] + x[1] * x[1];
+
+        // Convert domain block into hollow cylinder
+        if ( rsq > Rout * Rout || rsq < ( Rin - Wall_th ) * ( Rin - Wall_th ) )
+            return false;
+        // Leave remaining bottom wall particles and remove particles in between
+        // inner and outer cylinder
+        if ( x[2] > low_corner[2] + Wall_th && rsq > Rin * Rin &&
+             rsq < ( Rout - Wall_th ) * ( Rout - Wall_th ) )
+            return false;
+
+        return true;
+    };
+    // Container particles should be frozen, never updated.
+    auto particles = CabanaPD::createParticles<memory_space, model_type>(
+        exec_space(), low_corner, high_corner, num_cells, halo_width,
+        CabanaPD::BaseOutput{}, create_container, 0, true );
+
+    // Create powder.
+    auto create_powder = KOKKOS_LAMBDA( const int, const double x[3] )
+    {
+        double rsq = x[0] * x[0] + x[1] * x[1];
+
+        // Only create particles in between inner and outer cylinder.
+        if ( x[2] > low_corner[2] + Wall_th && rsq > ( Rin ) * ( Rin ) &&
+             rsq < ( Rout - Wall_th ) * ( Rout - Wall_th ) )
+            return true;
+
+        return false;
+    };
+    particles->createParticles( exec_space(), Cabana::InitRandom{},
+                                create_powder, particles->numFrozen() );
+
+    // Set density/volumes.
+    auto rho = particles->sliceDensity();
+    auto vol = particles->sliceVolume();
+    auto init_functor = KOKKOS_LAMBDA( const int pid )
+    {
+        rho( pid ) = rho0;
+        vol( pid ) = vol0;
+    };
+    particles->updateParticles( exec_space{}, init_functor );
+
+    // ====================================================
+    //                   Create solver
+    // ====================================================
+    auto cabana_pd = CabanaPD::createSolver<memory_space>( inputs, particles,
+                                                           contact_model );
+
+    // ====================================================
+    //                   Simulation init
+    // ====================================================
+    cabana_pd->init();
+
+    // Remove any points that are too close.
+    Kokkos::View<int*, memory_space> keep( "keep_points",
+                                           particles->numLocal() );
+    Kokkos::deep_copy( keep, 1 );
+    int num_keep;
+    auto num_frozen = particles->numFrozen();
+    auto f = particles->sliceForce();
+    auto remove_functor = KOKKOS_LAMBDA( const int pid, int& k )
+    {
+        auto f_mag = Kokkos::hypot( f( pid, 0 ), f( pid, 1 ), f( pid, 2 ) );
+        if ( f_mag > 1e8 )
+            keep( pid - num_frozen ) = 0;
+        else
+            k++;
+    };
+    Kokkos::RangePolicy<exec_space> policy( num_frozen,
+                                            particles->localOffset() );
+    Kokkos::parallel_reduce( "remove", policy, remove_functor,
+                             Kokkos::Sum<int>( num_keep ) );
+    cabana_pd->particles->remove( num_keep, keep );
+    // FIXME: Will need to rebuild ghosts.
+
+    // ====================================================
+    //                   Boundary condition
+    // ====================================================
+    f = cabana_pd->particles->sliceForce();
+    rho = cabana_pd->particles->sliceDensity();
+    auto body_functor = KOKKOS_LAMBDA( const int pid, const double )
+    {
+        f( pid, 2 ) -= 9.8 * rho( pid ); // * vol( pid );
+    };
+    auto gravity = CabanaPD::createBodyTerm( body_functor, true );
+
+    // ====================================================
+    //                   Simulation run
+    // ====================================================
+    cabana_pd->run( gravity );
+}
+
+// Initialize MPI+Kokkos.
+int main( int argc, char* argv[] )
+{
+    MPI_Init( &argc, &argv );
+    Kokkos::initialize( argc, argv );
+
+    powderSettlingExample( argv[1] );
+
+    Kokkos::finalize();
+    MPI_Finalize();
+}

examples/mechanics/fragmenting_cylinder.cpp

@@ -86,7 +86,8 @@ void fragmentingCylinderExample( const std::string filename )
     };
 
     auto particles = CabanaPD::createParticles<memory_space, model_type>(
-        exec_space(), low_corner, high_corner, num_cells, halo_width, init_op );
+        exec_space(), low_corner, high_corner, num_cells, Cabana::InitRandom{},
+        halo_width, init_op );
 
     auto rho = particles->sliceDensity();
     auto x = particles->sliceReferencePosition();
@@ -100,28 +101,12 @@ void fragmentingCylinderExample( const std::string filename )
 
     auto dx = particles->dx;
     double height = inputs["system_size"][2];
-    double factor = inputs["grid_perturbation_factor"];
 
-    using pool_type = Kokkos::Random_XorShift64_Pool<exec_space>;
-    using random_type = Kokkos::Random_XorShift64<exec_space>;
-    pool_type pool;
-    int seed = 456854;
-    pool.init( seed, particles->localOffset() );
     auto init_functor = KOKKOS_LAMBDA( const int pid )
     {
         // Density
         rho( pid ) = rho0;
 
-        // Perturb particle positions
-        auto gen = pool.get_state();
-        for ( std::size_t d = 0; d < 3; d++ )
-        {
-            auto rand =
-                Kokkos::rand<random_type, double>::draw( gen, 0.0, 1.0 );
-            x( pid, d ) += ( 2.0 * rand - 1.0 ) * factor * dx[d];
-        }
-        pool.free_state( gen );
-
         // Velocity
         double zfactor = ( ( x( pid, 2 ) - zmin ) / ( 0.5 * height ) ) - 1;
         double vr = vrmax - vrmin * zfactor * zfactor;