diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 27c4790a..bcbf987a 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,2 +1,3 @@ add_subdirectory(mechanics) add_subdirectory(thermomechanics) +add_subdirectory(dem) diff --git a/examples/dem/CMakeLists.txt b/examples/dem/CMakeLists.txt new file mode 100644 index 00000000..b2199147 --- /dev/null +++ b/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}) diff --git a/examples/dem/inputs/powder_fill.json b/examples/dem/inputs/powder_fill.json new file mode 100644 index 00000000..46136d9e --- /dev/null +++ b/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} +} diff --git a/examples/dem/powder_fill.cpp b/examples/dem/powder_fill.cpp new file mode 100644 index 00000000..b7ccd0db --- /dev/null +++ b/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 +#include + +#include "mpi.h" + +#include + +#include + +// 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 low_corner = inputs["low_corner"]; + std::array high_corner = inputs["high_corner"]; + std::array 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( + 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( inputs, particles, + contact_model ); + + // ==================================================== + // Simulation init + // ==================================================== + cabana_pd->init(); + + // Remove any points that are too close. + Kokkos::View 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 policy( num_frozen, + particles->localOffset() ); + Kokkos::parallel_reduce( "remove", policy, remove_functor, + Kokkos::Sum( 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(); +} diff --git a/examples/mechanics/fragmenting_cylinder.cpp b/examples/mechanics/fragmenting_cylinder.cpp index 12e86cd6..60012124 100644 --- a/examples/mechanics/fragmenting_cylinder.cpp +++ b/examples/mechanics/fragmenting_cylinder.cpp @@ -86,7 +86,8 @@ void fragmentingCylinderExample( const std::string filename ) }; auto particles = CabanaPD::createParticles( - 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; - using random_type = Kokkos::Random_XorShift64; - 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::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; diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 08d69a7b..71148753 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -173,7 +173,8 @@ class Particles , _plist_f( "forces" ) { createDomain( low_corner, high_corner, num_cells ); - createParticles( exec_space, num_previous, create_frozen ); + createParticles( exec_space, Cabana::InitUniform{}, num_previous, + create_frozen ); } // Constructor which initializes particles on regular grid with @@ -189,7 +190,42 @@ class Particles , _plist_f( "forces" ) { createDomain( low_corner, high_corner, num_cells ); - createParticles( exec_space, user_create, num_previous, create_frozen ); + createParticles( exec_space, Cabana::InitUniform{}, user_create, + num_previous, create_frozen ); + } + + // Constructor which initializes particles on regular grid, randomly per + // cell. + template + Particles( const ExecSpace& exec_space, std::array low_corner, + std::array high_corner, + const std::array num_cells, const int max_halo_width, + Cabana::InitRandom random, const std::size_t num_previous = 0, + const bool create_frozen = false ) + : halo_width( max_halo_width ) + , _plist_x( "positions" ) + , _plist_f( "forces" ) + { + createDomain( low_corner, high_corner, num_cells ); + createParticles( exec_space, random, num_previous, create_frozen ); + } + + // Constructor which initializes particles on regular grid with + // customization, randomly per cell. + template + Particles( const ExecSpace& exec_space, std::array low_corner, + std::array high_corner, + const std::array num_cells, const int max_halo_width, + Cabana::InitRandom random, UserFunctor user_create, + const std::size_t num_previous = 0, + const bool create_frozen = false ) + : halo_width( max_halo_width ) + , _plist_x( "positions" ) + , _plist_f( "forces" ) + { + createDomain( low_corner, high_corner, num_cells ); + createParticles( exec_space, random, user_create, num_previous, + create_frozen ); } // Constructor with existing particle data. @@ -252,22 +288,34 @@ class Particles _init_timer.stop(); } - template - void createParticles( const ExecSpace& exec_space, - const std::size_t num_previous = 0, - const bool create_frozen = false ) + template + void + createParticles( const ExecSpace& exec_space, InitType init_type, + const std::size_t num_previous = 0, + const bool create_frozen = false, + typename std::enable_if< + ( std::is_same::value || + std::is_same::value ), + int>::type* = 0 ) { auto empty = KOKKOS_LAMBDA( const int, const double[dim] ) { return true; }; - createParticles( exec_space, empty, num_previous, create_frozen ); + createParticles( exec_space, init_type, empty, num_previous, + create_frozen ); } - template - void createParticles( const ExecSpace& exec_space, UserFunctor user_create, - const std::size_t num_previous = 0, - const bool create_frozen = false ) + template + void + createParticles( const ExecSpace& exec_space, InitType init_type, + UserFunctor user_create, + const std::size_t num_previous = 0, + const bool create_frozen = false, + typename std::enable_if< + ( std::is_same::value || + std::is_same::value ), + int>::type* = 0 ) { _init_timer.start(); // Create a local mesh and owned space. @@ -320,15 +368,11 @@ class Particles return create; }; - local_offset = Cabana::Grid::createParticles( - Cabana::InitUniform{}, exec_space, create_functor, _plist_x, 1, - *local_grid, num_previous ); - resize( local_offset, 0 ); - size = _plist_x.size(); - - // Only set this value if this generation of particles should be frozen. - if ( create_frozen ) - frozen_offset = size; + // Fence inside create. + auto local_created = Cabana::Grid::createParticles( + init_type, exec_space, create_functor, _plist_x, particles_per_cell, + *local_grid, num_previous, false ); + resize( local_created, 0, create_frozen ); updateGlobal(); _init_timer.stop(); @@ -336,18 +380,20 @@ class Particles // Store custom created particle positions and volumes. template - void createParticles( const ExecSpace, const PositionType& x, - const VolumeType& vol, - const std::size_t num_previous = 0, - const bool create_frozen = false ) + void createParticles( + const ExecSpace, const PositionType& x, const VolumeType& vol, + const std::size_t num_previous = 0, const bool create_frozen = false, + typename std::enable_if<( Cabana::is_slice::value || + Kokkos::is_view::value ) && + ( Cabana::is_slice::value || + Kokkos::is_view::value ), + int>::type* = 0 ) { // Ensure valid previous particles. assert( num_previous <= referenceOffset() ); // Ensure matching input sizes. assert( vol.size() == x.extent( 0 ) ); - resize( vol.size() + num_previous, 0 ); - if ( create_frozen ) - frozen_offset = size; + resize( vol.size() + num_previous, 0, create_frozen ); auto p_x = sliceReferencePosition(); auto p_vol = sliceVolume(); @@ -381,6 +427,7 @@ class Particles nofail( pid ) = 0; rho( pid ) = 1.0; } ); + Kokkos::fence(); updateGlobal(); } @@ -406,6 +453,7 @@ class Particles Kokkos::parallel_for( "CabanaPD::Particles::update_particles", policy, KOKKOS_LAMBDA( const int pid ) { init_functor( pid ); } ); + Kokkos::fence(); _timer.stop(); } @@ -510,9 +558,13 @@ class Particles //_timer.stop(); } - void resize( int new_local, int new_ghost ) + void resize( int new_local, int new_ghost, + const bool create_frozen = false ) { _timer.start(); + if ( new_ghost > 0 ) + assert( create_frozen == false ); + local_offset = new_local; num_ghost = new_ghost; size = new_local + new_ghost; @@ -524,10 +576,46 @@ class Particles _plist_f.aosoa().resize( localOffset() ); _aosoa_other.resize( localOffset() ); _aosoa_nofail.resize( referenceOffset() ); + size = _plist_x.size(); + if ( create_frozen ) + frozen_offset = size; + _timer.stop(); + }; + + void shrink() + { + _timer.start(); + _plist_x.aosoa().shrinkToFit(); + _aosoa_u.shrinkToFit(); + _aosoa_y.shrinkToFit(); + _aosoa_vol.shrinkToFit(); + _plist_f.aosoa().shrinkToFit(); + _aosoa_other.shrinkToFit(); + _aosoa_nofail.shrinkToFit(); _timer.stop(); }; + template + void remove( const int num_keep, const KeepType& keep ) + { + Cabana::remove( execution_space(), num_keep, keep, _plist_x.aosoa(), + numFrozen() ); + Cabana::remove( execution_space(), num_keep, keep, _plist_f.aosoa(), + numFrozen() ); + Cabana::remove( execution_space(), num_keep, keep, _aosoa_u, + numFrozen() ); + Cabana::remove( execution_space(), num_keep, keep, _aosoa_vol, + numFrozen() ); + Cabana::remove( execution_space(), num_keep, keep, _aosoa_y, + numFrozen() ); + Cabana::remove( execution_space(), num_keep, keep, _aosoa_other, + numFrozen() ); + Cabana::remove( execution_space(), num_keep, keep, _aosoa_nofail, + numFrozen() ); + resize( frozen_offset + num_keep, 0 ); + } + auto getPosition( const bool use_reference ) { if ( use_reference ) @@ -643,13 +731,9 @@ class Particles _init_timer.stop(); } - // Constructor which initializes particles on regular grid. - template - Particles( const ExecSpace& exec_space, std::array low_corner, - std::array high_corner, - const std::array num_cells, const int max_halo_width ) - : base_type( exec_space, low_corner, high_corner, num_cells, - max_halo_width ) + template + Particles( Args&&... args ) + : base_type( std::forward( args )... ) { _init_timer.start(); _aosoa_m = aosoa_m_type( "Particle Weighted Volumes", @@ -771,12 +855,9 @@ class Particles } // Constructor which initializes particles on regular grid. - template - Particles( const ExecSpace& exec_space, std::array low_corner, - std::array high_corner, - const std::array num_cells, const int max_halo_width ) - : base_type( exec_space, low_corner, high_corner, num_cells, - max_halo_width ) + template + Particles( Args&&... args ) + : base_type( std::forward( args )... ) { _aosoa_temp = aosoa_temp_type( "Particle Temperature", base_type::localOffset() ); @@ -816,9 +897,10 @@ class Particles return temp_a; } - void resize( int new_local, int new_ghost ) + template + void resize( Args&&... args ) { - base_type::resize( new_local, new_ghost ); + base_type::resize( std::forward( args )... ); _aosoa_temp.resize( base_type::referenceOffset() ); } @@ -922,9 +1004,10 @@ class Particles return Cabana::slice<1>( _aosoa_output, "damage" ); } - void resize( int new_local, int new_ghost ) + template + void resize( Args&&... args ) { - base_type::resize( new_local, new_ghost ); + base_type::resize( std::forward( args )... ); _aosoa_output.resize( base_type::localOffset() ); } @@ -993,14 +1076,16 @@ template low_corner, std::array high_corner, const std::array num_cells, - const int max_halo_width, OutputType, + const int max_halo_width, OutputType, const std::size_t num_previous = 0, + const bool create_frozen = false, typename std::enable_if<( is_temperature_dependent::value && is_output::value ), int>::type* = 0 ) { return std::make_shared>( - exec_space, low_corner, high_corner, num_cells, max_halo_width ); + exec_space, low_corner, high_corner, num_cells, max_halo_width, + num_previous, create_frozen ); } // Backwards compatible versions with energy output by default. @@ -1030,12 +1115,14 @@ auto createParticles( const ExecSpace& exec_space, std::array low_corner, std::array high_corner, const std::array num_cells, const int max_halo_width, UserFunctor user, OutputType, + const std::size_t num_previous = 0, const bool create_frozen = false, typename std::enable_if<( is_output::value ), int>::type* = 0 ) { return std::make_shared< CabanaPD::Particles>( - exec_space, low_corner, high_corner, num_cells, max_halo_width, user ); + exec_space, low_corner, high_corner, num_cells, max_halo_width, user, + num_previous, create_frozen ); } template low_corner, std::array high_corner, const std::array num_cells, - const int max_halo_width, UserFunctor user ) + Cabana::InitRandom random, const int max_halo_width, + UserFunctor user, const std::size_t num_previous = 0, + const bool create_frozen = false ) { - return createParticles( - exec_space, low_corner, high_corner, num_cells, max_halo_width, user, - EnergyOutput{} ); + return std::make_shared< + CabanaPD::Particles>( + exec_space, low_corner, high_corner, num_cells, max_halo_width, random, + user, num_previous, create_frozen ); } template +auto createParticles( const ExecSpace& exec_space, + std::array low_corner, + std::array high_corner, + const std::array num_cells, + const int max_halo_width, UserFunctor user_create, + const std::size_t num_previous = 0, + const bool create_frozen = false ) +{ + return std::make_shared< + CabanaPD::Particles>( + exec_space, low_corner, high_corner, num_cells, max_halo_width, + user_create, num_previous, create_frozen ); +} + +template +auto createParticles( + const ExecSpace& exec_space, std::array low_corner, + std::array high_corner, const std::array num_cells, + const int max_halo_width, OutputType, UserFunctor user_create, + const std::size_t num_previous = 0, const bool create_frozen = false ) +{ + return std::make_shared< + CabanaPD::Particles>( + exec_space, low_corner, high_corner, num_cells, max_halo_width, + user_create, num_previous, create_frozen ); +} + template diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 2696605f..b5e9ab79 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -109,8 +109,8 @@ class Solver Solver( input_type _inputs, std::shared_ptr _particles, force_model_type force_model ) - : inputs( _inputs ) - , particles( _particles ) + : particles( _particles ) + , inputs( _inputs ) , _init_time( 0.0 ) { setup( force_model ); @@ -118,8 +118,8 @@ class Solver Solver( input_type _inputs, std::shared_ptr _particles, force_model_type force_model, contact_model_type contact_model ) - : inputs( _inputs ) - , particles( _particles ) + : particles( _particles ) + , inputs( _inputs ) , _init_time( 0.0 ) { setup( force_model ); @@ -473,6 +473,8 @@ class Solver bool output_reference; double dt; int thermal_subcycle_steps; + // Sometimes necessary to update particles after solver creation. + std::shared_ptr particles; protected: template @@ -489,7 +491,6 @@ class Solver // Core modules. input_type inputs; - std::shared_ptr particles; std::shared_ptr comm; std::shared_ptr integrator; std::shared_ptr force; diff --git a/unit_test/tstParticles.hpp b/unit_test/tstParticles.hpp index d188cdcc..d7204c32 100644 --- a/unit_test/tstParticles.hpp +++ b/unit_test/tstParticles.hpp @@ -112,7 +112,7 @@ void testCreateFrozenParticles() return false; }; // Create more, starting from the current number of frozen points. - particles.createParticles( exec_space{}, init_top, + particles.createParticles( exec_space{}, Cabana::InitUniform{}, init_top, particles.frozenOffset() ); // Check expected values for each block of particles.