From d1dd750d91e8d9895214cad6754655b5c58376e5 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 30 Dec 2024 10:29:45 -0500 Subject: [PATCH 01/11] Initial powder fill example --- examples/CMakeLists.txt | 1 + examples/dem/CMakeLists.txt | 3 + examples/dem/inputs/powder_fill.json | 19 ++++ examples/dem/powder_fill.cpp | 146 +++++++++++++++++++++++++++ 4 files changed, 169 insertions(+) create mode 100644 examples/dem/CMakeLists.txt create mode 100644 examples/dem/inputs/powder_fill.json create mode 100644 examples/dem/powder_fill.cpp 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..196b2bdc --- /dev/null +++ b/examples/dem/inputs/powder_fill.json @@ -0,0 +1,19 @@ +{ + "num_cells" : {"value": [201, 201, 257]}, + "system_size" : {"value": [0.23749, 0.23749, 0.3048], "unit": "m"}, + "density" : {"value": 7.95e3, "unit": "kg/m^3"}, + "volume" : {"value": 5.236e-13, "unit": "m^3"}, + "elastic_modulus" : {"value": 195.6e9, "unit": "Pa"}, + "poisson_ratio" : {"value": 0.25, "unit": ""}, + "restitution" : {"value": 0.5}, + "radius" : {"value": 5e-5 }, + "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": 4e-03, "unit": "s"}, + "timestep" : {"value": 4e-07, "unit": "s"}, + "timestep_safety_factor" : {"value": 0.9}, + "output_frequency" : {"value": 100}, + "output_reference" : {"value": false} +} diff --git a/examples/dem/powder_fill.cpp b/examples/dem/powder_fill.cpp new file mode 100644 index 00000000..616c0429 --- /dev/null +++ b/examples/dem/powder_fill.cpp @@ -0,0 +1,146 @@ +/**************************************************************************** + * 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, 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, + create_container, 0, true ); + std::cout << particles->numFrozen() << " " << particles->numLocal() + << std::endl; + + // Create powder. + auto dx = particles->dx[0] * 2.0; + 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 + dx && + rsq > ( Rin + dx ) * ( Rin + dx ) && + rsq < ( Rout - Wall_th - dx ) * ( Rout - Wall_th - dx ) ) + return true; + + return false; + }; + particles->createParticles( exec_space(), Cabana::InitRandom{}, + create_powder, particles->numFrozen() ); + + auto rho = particles->sliceDensity(); + auto vol = particles->sliceVolume(); + auto init_functor = KOKKOS_LAMBDA( const int pid ) + { + // Density + rho( pid ) = rho0; + vol( pid ) = vol0; + }; + particles->updateParticles( exec_space{}, init_functor ); + + // ==================================================== + // Create solver + // ==================================================== + auto cabana_pd = CabanaPD::createSolver( inputs, particles, + contact_model ); + + // ==================================================== + // Simulation run + // ==================================================== + cabana_pd->init(); + cabana_pd->run(); +} + +// Initialize MPI+Kokkos. +int main( int argc, char* argv[] ) +{ + MPI_Init( &argc, &argv ); + Kokkos::initialize( argc, argv ); + + powderSettlingExample( argv[1] ); + + Kokkos::finalize(); + MPI_Finalize(); +} From 13a4fcfa85e4dc80c397dc5aaf1dedcc593b768b Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 30 Dec 2024 10:40:35 -0500 Subject: [PATCH 02/11] Enable uniform or random particle init --- examples/dem/powder_fill.cpp | 2 +- examples/mechanics/fragmenting_cylinder.cpp | 19 +-- src/CabanaPD_Particles.hpp | 170 ++++++++++++++------ unit_test/tstParticles.hpp | 2 +- 4 files changed, 124 insertions(+), 69 deletions(-) diff --git a/examples/dem/powder_fill.cpp b/examples/dem/powder_fill.cpp index 616c0429..abe4fc51 100644 --- a/examples/dem/powder_fill.cpp +++ b/examples/dem/powder_fill.cpp @@ -89,7 +89,7 @@ void powderSettlingExample( const std::string filename ) // Container particles should be frozen, never updated. auto particles = CabanaPD::createParticles( exec_space(), low_corner, high_corner, num_cells, halo_width, - create_container, 0, true ); + create_container, CabanaPD::BaseOutput{}, 0, true ); std::cout << particles->numFrozen() << " " << particles->numLocal() << std::endl; 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..7cd7c9b5 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,10 @@ 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; + 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 +379,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(); @@ -510,9 +555,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,7 +573,10 @@ 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(); }; @@ -643,13 +695,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 +819,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 +861,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 +968,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 +1040,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 +1079,31 @@ 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 +auto createParticles( const ExecSpace& exec_space, + std::array low_corner, + std::array high_corner, + const std::array num_cells, + Cabana::InitRandom random, const int max_halo_width, + UserFunctor user, 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, random, + user, num_previous, create_frozen ); } template low_corner, std::array high_corner, const std::array num_cells, - const int max_halo_width, UserFunctor user ) + 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{} ); + EnergyOutput{}, num_previous, create_frozen ); } template Date: Mon, 30 Dec 2024 10:42:22 -0500 Subject: [PATCH 03/11] Add options to remove particles or shrink memory --- src/CabanaPD_Particles.hpp | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 7cd7c9b5..00694d38 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -580,6 +580,39 @@ class Particles _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 ) From 6014a55b2e5706e005499854cbd0f8fbb4f7b904 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 30 Dec 2024 10:42:40 -0500 Subject: [PATCH 04/11] TODO: Add all Particles creation methods --- src/CabanaPD_Particles.hpp | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 00694d38..507adfea 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -1168,6 +1168,23 @@ auto createParticles( EnergyOutput{} ); } +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 From 1350aab7ba5b232acfd8efc6c44cafedf111f81d Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 30 Dec 2024 14:25:17 -0500 Subject: [PATCH 05/11] Remove particles that are too close --- examples/dem/powder_fill.cpp | 34 +++++++++++++++++++++++++++++++++- 1 file changed, 33 insertions(+), 1 deletion(-) diff --git a/examples/dem/powder_fill.cpp b/examples/dem/powder_fill.cpp index abe4fc51..8aac2b4b 100644 --- a/examples/dem/powder_fill.cpp +++ b/examples/dem/powder_fill.cpp @@ -127,9 +127,41 @@ void powderSettlingExample( const std::string filename ) contact_model ); // ==================================================== - // Simulation run + // 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( pid, 0 ) > 0 ) + std::cout << f_mag << std::endl; + 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 ) ); + std::cout << num_frozen << " " << particles->localOffset() << " " + << particles->numLocal() - num_keep << " " << keep.size() + << std::endl; + cabana_pd->particles->remove( num_keep, keep ); + std::cout << particles->numFrozen() << " " << particles->numLocal() + << std::endl; + // Need to rebuild ghosts.. + // ==================================================== + // Simulation run + // ==================================================== cabana_pd->run(); } From a270efec19cfcf3958d444fb016162cb4500582c Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 30 Dec 2024 14:25:58 -0500 Subject: [PATCH 06/11] Misc changes to run powder settling --- examples/dem/inputs/powder_fill.json | 12 +++++------ examples/dem/powder_fill.cpp | 26 +++++++++++++++++------- src/CabanaPD_Particles.hpp | 30 ++++++++++++++-------------- src/CabanaPD_Solver.hpp | 3 ++- 4 files changed, 42 insertions(+), 29 deletions(-) diff --git a/examples/dem/inputs/powder_fill.json b/examples/dem/inputs/powder_fill.json index 196b2bdc..672c0c9a 100644 --- a/examples/dem/inputs/powder_fill.json +++ b/examples/dem/inputs/powder_fill.json @@ -1,9 +1,9 @@ { - "num_cells" : {"value": [201, 201, 257]}, - "system_size" : {"value": [0.23749, 0.23749, 0.3048], "unit": "m"}, + "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-13, "unit": "m^3"}, - "elastic_modulus" : {"value": 195.6e9, "unit": "Pa"}, + "elastic_modulus" : {"value": 195.6e6, "unit": "Pa"}, "poisson_ratio" : {"value": 0.25, "unit": ""}, "restitution" : {"value": 0.5}, "radius" : {"value": 5e-5 }, @@ -11,9 +11,9 @@ "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": 4e-03, "unit": "s"}, - "timestep" : {"value": 4e-07, "unit": "s"}, + "final_time" : {"value": 4e-01, "unit": "s"}, + "timestep" : {"value": 2e-07, "unit": "s"}, "timestep_safety_factor" : {"value": 0.9}, - "output_frequency" : {"value": 100}, + "output_frequency" : {"value": 1}, "output_reference" : {"value": false} } diff --git a/examples/dem/powder_fill.cpp b/examples/dem/powder_fill.cpp index 8aac2b4b..d6e8badb 100644 --- a/examples/dem/powder_fill.cpp +++ b/examples/dem/powder_fill.cpp @@ -59,7 +59,7 @@ void powderSettlingExample( const std::string filename ) // Force model // ==================================================== using model_type = CabanaPD::HertzianModel; - model_type contact_model( delta, radius, nu, E, e ); + model_type contact_model( delta / m * 0.9, radius, nu, E, e ); // ==================================================== // Custom particle initialization @@ -89,7 +89,7 @@ void powderSettlingExample( const std::string filename ) // Container particles should be frozen, never updated. auto particles = CabanaPD::createParticles( exec_space(), low_corner, high_corner, num_cells, halo_width, - create_container, CabanaPD::BaseOutput{}, 0, true ); + CabanaPD::BaseOutput{}, create_container, 0, true ); std::cout << particles->numFrozen() << " " << particles->numLocal() << std::endl; @@ -100,21 +100,22 @@ void powderSettlingExample( const std::string filename ) 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 + dx && - rsq > ( Rin + dx ) * ( Rin + dx ) && - rsq < ( Rout - Wall_th - dx ) * ( Rout - Wall_th - dx ) ) + 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() ); + std::cout << particles->numFrozen() << " " << particles->numLocal() + << std::endl; + // Set density/volumes. auto rho = particles->sliceDensity(); auto vol = particles->sliceVolume(); auto init_functor = KOKKOS_LAMBDA( const int pid ) { - // Density rho( pid ) = rho0; vol( pid ) = vol0; }; @@ -159,10 +160,21 @@ void powderSettlingExample( const std::string filename ) std::cout << particles->numFrozen() << " " << particles->numLocal() << std::endl; // Need to rebuild ghosts.. + + // ==================================================== + // Boundary condition + // ==================================================== + f = cabana_pd->particles->sliceForce(); + 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(); + cabana_pd->run( gravity ); } // Initialize MPI+Kokkos. diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 507adfea..198cffe2 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -1139,21 +1139,6 @@ auto createParticles( const ExecSpace& exec_space, 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, - 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{}, num_previous, create_frozen ); -} - template auto createParticles( @@ -1185,6 +1170,21 @@ auto createParticles( const ExecSpace& exec_space, 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..461d5905 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -474,6 +474,8 @@ class Solver double dt; int thermal_subcycle_steps; + std::shared_ptr particles; + protected: template void init_prenotch( Prenotch prenotch ) @@ -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; From 387a8ca157d6e5b0bff81823465e73b8182361d1 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 30 Dec 2024 16:25:45 -0500 Subject: [PATCH 07/11] fixup: warnings --- examples/dem/powder_fill.cpp | 13 +------------ src/CabanaPD_Solver.hpp | 10 +++++----- 2 files changed, 6 insertions(+), 17 deletions(-) diff --git a/examples/dem/powder_fill.cpp b/examples/dem/powder_fill.cpp index d6e8badb..4de164d4 100644 --- a/examples/dem/powder_fill.cpp +++ b/examples/dem/powder_fill.cpp @@ -90,8 +90,6 @@ void powderSettlingExample( const std::string filename ) auto particles = CabanaPD::createParticles( exec_space(), low_corner, high_corner, num_cells, halo_width, CabanaPD::BaseOutput{}, create_container, 0, true ); - std::cout << particles->numFrozen() << " " << particles->numLocal() - << std::endl; // Create powder. auto dx = particles->dx[0] * 2.0; @@ -108,8 +106,6 @@ void powderSettlingExample( const std::string filename ) }; particles->createParticles( exec_space(), Cabana::InitRandom{}, create_powder, particles->numFrozen() ); - std::cout << particles->numFrozen() << " " << particles->numLocal() - << std::endl; // Set density/volumes. auto rho = particles->sliceDensity(); @@ -142,8 +138,6 @@ void powderSettlingExample( const std::string filename ) 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( pid, 0 ) > 0 ) - std::cout << f_mag << std::endl; if ( f_mag > 1e8 ) keep( pid - num_frozen ) = 0; else @@ -153,13 +147,8 @@ void powderSettlingExample( const std::string filename ) particles->localOffset() ); Kokkos::parallel_reduce( "remove", policy, remove_functor, Kokkos::Sum( num_keep ) ); - std::cout << num_frozen << " " << particles->localOffset() << " " - << particles->numLocal() - num_keep << " " << keep.size() - << std::endl; cabana_pd->particles->remove( num_keep, keep ); - std::cout << particles->numFrozen() << " " << particles->numLocal() - << std::endl; - // Need to rebuild ghosts.. + // FIXME: Will need to rebuild ghosts. // ==================================================== // Boundary condition diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 461d5905..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,7 +473,7 @@ class Solver bool output_reference; double dt; int thermal_subcycle_steps; - + // Sometimes necessary to update particles after solver creation. std::shared_ptr particles; protected: From 65c106729c43b49337d5715f2c47ddf6987d82a0 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 30 Dec 2024 16:29:53 -0500 Subject: [PATCH 08/11] fixup: unused --- examples/dem/powder_fill.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/dem/powder_fill.cpp b/examples/dem/powder_fill.cpp index 4de164d4..d1c78f8d 100644 --- a/examples/dem/powder_fill.cpp +++ b/examples/dem/powder_fill.cpp @@ -92,7 +92,6 @@ void powderSettlingExample( const std::string filename ) CabanaPD::BaseOutput{}, create_container, 0, true ); // Create powder. - auto dx = particles->dx[0] * 2.0; auto create_powder = KOKKOS_LAMBDA( const int, const double x[3] ) { double rsq = x[0] * x[0] + x[1] * x[1]; From fb2ba6ffc0926180254dd77b4a8425253fac00e5 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 3 Jan 2025 10:24:58 -0500 Subject: [PATCH 09/11] fixup: fences --- src/CabanaPD_Particles.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 198cffe2..71148753 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -368,6 +368,7 @@ class Particles return create; }; + // 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 ); @@ -426,6 +427,7 @@ class Particles nofail( pid ) = 0; rho( pid ) = 1.0; } ); + Kokkos::fence(); updateGlobal(); } @@ -451,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(); } From a8ce7be96c5f19d4742bc342ad1197de2b954553 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 15 Jan 2025 10:28:15 -0500 Subject: [PATCH 10/11] fixup: missing slice update --- examples/dem/powder_fill.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/dem/powder_fill.cpp b/examples/dem/powder_fill.cpp index d1c78f8d..b7ccd0db 100644 --- a/examples/dem/powder_fill.cpp +++ b/examples/dem/powder_fill.cpp @@ -153,6 +153,7 @@ void powderSettlingExample( const std::string filename ) // 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 ); From 44dfaa27ab2b21b8c8e9e60a05cc2c4aeca99a58 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Thu, 16 Jan 2025 17:23:57 -0500 Subject: [PATCH 11/11] fixup: Austin's changes --- examples/dem/inputs/powder_fill.json | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/dem/inputs/powder_fill.json b/examples/dem/inputs/powder_fill.json index 672c0c9a..46136d9e 100644 --- a/examples/dem/inputs/powder_fill.json +++ b/examples/dem/inputs/powder_fill.json @@ -2,18 +2,18 @@ "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-13, "unit": "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-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": 4e-01, "unit": "s"}, - "timestep" : {"value": 2e-07, "unit": "s"}, + "final_time" : {"value": 4e0, "unit": "s"}, + "timestep" : {"value": 2e-6, "unit": "s"}, "timestep_safety_factor" : {"value": 0.9}, - "output_frequency" : {"value": 1}, + "output_frequency" : {"value": 1000}, "output_reference" : {"value": false} }