From 27a100b3810936b8e7d83311d0e3b0dfaba8ff29 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 18 Nov 2024 23:23:54 -0500 Subject: [PATCH 01/23] Add PD contact and impactor example --- examples/mechanics/CMakeLists.txt | 5 +- examples/mechanics/crack_branching.cpp | 4 +- examples/mechanics/disk_impact.cpp | 101 ++++++++++++ examples/mechanics/inputs/disk_impact.json | 16 ++ examples/mechanics/kalthoff_winkler.cpp | 4 +- src/CabanaPD.hpp | 1 + src/CabanaPD_Contact.hpp | 174 +++++++++++++++++++++ src/CabanaPD_Force.hpp | 15 +- src/CabanaPD_Solver.hpp | 130 +++++++++++---- src/CabanaPD_Types.hpp | 23 ++- 10 files changed, 431 insertions(+), 42 deletions(-) create mode 100644 examples/mechanics/disk_impact.cpp create mode 100644 examples/mechanics/inputs/disk_impact.json create mode 100644 src/CabanaPD_Contact.hpp diff --git a/examples/mechanics/CMakeLists.txt b/examples/mechanics/CMakeLists.txt index a10b05f1..19148fa3 100644 --- a/examples/mechanics/CMakeLists.txt +++ b/examples/mechanics/CMakeLists.txt @@ -10,4 +10,7 @@ target_link_libraries(CrackBranching LINK_PUBLIC CabanaPD) add_executable(FragmentingCylinder fragmenting_cylinder.cpp) target_link_libraries(FragmentingCylinder LINK_PUBLIC CabanaPD) -install(TARGETS ElasticWave KalthoffWinkler CrackBranching FragmentingCylinder DESTINATION ${CMAKE_INSTALL_BINDIR}) +add_executable(DiskImpact disk_impact.cpp) +target_link_libraries(DiskImpact LINK_PUBLIC CabanaPD) + +install(TARGETS ElasticWave KalthoffWinkler CrackBranching FragmentingCylinder DiskImpact DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/examples/mechanics/crack_branching.cpp b/examples/mechanics/crack_branching.cpp index 65a8edb8..b380472a 100644 --- a/examples/mechanics/crack_branching.cpp +++ b/examples/mechanics/crack_branching.cpp @@ -114,7 +114,7 @@ void crackBranchingExample( const std::string filename ) // Create solver // ==================================================== auto cabana_pd = CabanaPD::createSolverFracture( - inputs, particles, force_model, prenotch ); + inputs, particles, force_model ); // ==================================================== // Boundary conditions @@ -137,7 +137,7 @@ void crackBranchingExample( const std::string filename ) // ==================================================== // Simulation run // ==================================================== - cabana_pd->init(); + cabana_pd->init( bc, prenotch ); cabana_pd->run( bc ); } diff --git a/examples/mechanics/disk_impact.cpp b/examples/mechanics/disk_impact.cpp new file mode 100644 index 00000000..0877b4a2 --- /dev/null +++ b/examples/mechanics/disk_impact.cpp @@ -0,0 +1,101 @@ +#include +#include +#include + +#include "mpi.h" + +#include + +#include + +// Simulate a sphere impacting a thin cylindrical plate. +void diskImpactExample( const std::string filename ) +{ + using exec_space = Kokkos::DefaultExecutionSpace; + using memory_space = typename exec_space::memory_space; + + CabanaPD::Inputs inputs( filename ); + + double K = inputs["bulk_modulus"]; + double G0 = inputs["fracture_energy"]; + double delta = inputs["horizon"]; + delta += 1e-10; + // Choose force model type. + using model_type = CabanaPD::ForceModel; + model_type force_model( delta, K, G0 ); + + // Create particles from mesh. + auto particles = CabanaPD::createParticles( + exec_space(), inputs ); + + std::array low_corner = inputs["low_corner"]; + std::array high_corner = inputs["high_corner"]; + double global_mesh_ext = high_corner[0] - low_corner[0]; + auto x = particles->sliceReferencePosition(); + auto create_func = KOKKOS_LAMBDA( const int, const double px[3] ) + { + auto width = global_mesh_ext / 2.0; + auto r2 = px[0] * px[0] + px[1] * px[1]; + if ( r2 > width * width ) + return false; + return true; + }; + particles->createParticles( exec_space{}, create_func ); + + // Define particle initialization. + auto v = particles->sliceVelocity(); + auto f = particles->sliceForce(); + auto rho = particles->sliceDensity(); + + double rho0 = inputs["density"]; + 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 ); + + double dx = particles->dx[0]; + double r_c = dx * 2.0; + CabanaPD::NormalRepulsionModel contact_model( delta, r_c, K ); + + // FIXME: use createSolver to switch backend at runtime. + auto cabana_pd = CabanaPD::createSolverFracture( + inputs, particles, force_model, contact_model ); + + double impact_r = inputs["impactor_radius"]; + double impact_v = inputs["impactor_velocity"]; + double impact_x = 0.0; + double impact_y = 0.0; + auto body_func = KOKKOS_LAMBDA( const int p, const double t ) + { + auto z = t * impact_v + impact_r; + double r = sqrt( ( x( p, 0 ) - impact_x ) * ( x( p, 0 ) - impact_x ) + + ( x( p, 1 ) - impact_y ) * ( x( p, 1 ) - impact_y ) + + ( x( p, 2 ) - z ) * ( x( p, 2 ) - z ) ); + if ( r < impact_r ) + { + double fmag = -1.0e17 * ( r - impact_r ) * ( r - impact_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; + } + }; + auto body = CabanaPD::createBodyTerm( body_func, true ); + + cabana_pd->init( body ); + cabana_pd->run( body ); +} + +int main( int argc, char* argv[] ) +{ + MPI_Init( &argc, &argv ); + + Kokkos::initialize( argc, argv ); + + diskImpactExample( argv[1] ); + + Kokkos::finalize(); + MPI_Finalize(); +} diff --git a/examples/mechanics/inputs/disk_impact.json b/examples/mechanics/inputs/disk_impact.json new file mode 100644 index 00000000..6a273f0b --- /dev/null +++ b/examples/mechanics/inputs/disk_impact.json @@ -0,0 +1,16 @@ +{ + "num_cells" : {"value": [149, 149, 5]}, + "low_corner" : {"value": [-0.037, -0.037, 2.5e-3], "unit": "m"}, + "high_corner" : {"value": [0.037, 0.037, 0.0], "unit": "m"}, + "density" : {"value": 2200.0, "unit": "kg/m^3"}, + "bulk_modulus" : {"value": 14.9e+9, "unit": "Pa"}, + "fracture_energy" : {"value": 10.0575, "unit": "J/m^2"}, + "horizon" : {"value": 0.0015, "unit": "m"}, + "final_time" : {"value": 2.0e-4, "unit": "s"}, + "timestep" : {"value": 1.0e-7, "unit": "s"}, + "timestep_safety_factor" : {"value": 0.85}, + "output_frequency" : {"value": 50}, + "output_reference" : {"value": false}, + "impactor_radius" : {"value": 5e-3, "unit": "m"}, + "impactor_velocity" : {"value": -100.0, "unit": "m/s"} +} diff --git a/examples/mechanics/kalthoff_winkler.cpp b/examples/mechanics/kalthoff_winkler.cpp index 78dab593..b6a3372b 100644 --- a/examples/mechanics/kalthoff_winkler.cpp +++ b/examples/mechanics/kalthoff_winkler.cpp @@ -116,7 +116,7 @@ void kalthoffWinklerExample( const std::string filename ) // Create solver // ==================================================== auto cabana_pd = CabanaPD::createSolverFracture( - inputs, particles, force_model, prenotch ); + inputs, particles, force_model ); // ==================================================== // Boundary conditions @@ -132,7 +132,7 @@ void kalthoffWinklerExample( const std::string filename ) // ==================================================== // Simulation run // ==================================================== - cabana_pd->init(); + cabana_pd->init( bc, prenotch ); cabana_pd->run( bc ); } diff --git a/src/CabanaPD.hpp b/src/CabanaPD.hpp index 7da7ad22..380ab056 100644 --- a/src/CabanaPD.hpp +++ b/src/CabanaPD.hpp @@ -16,6 +16,7 @@ #include #include #include +#include #include #include #include diff --git a/src/CabanaPD_Contact.hpp b/src/CabanaPD_Contact.hpp new file mode 100644 index 00000000..40cd81f8 --- /dev/null +++ b/src/CabanaPD_Contact.hpp @@ -0,0 +1,174 @@ +/**************************************************************************** + * 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 * + ****************************************************************************/ + +#ifndef CONTACT_H +#define CONTACT_H + +#include + +#include +#include +#include + +namespace CabanaPD +{ +/****************************************************************************** + Contact model +******************************************************************************/ +struct ContactModel +{ + double delta; + double Rc; + + ContactModel(){}; + // PD horizon + // Contact radius + ContactModel( const double _delta, const double _Rc ) + : delta( _delta ) + , Rc( _Rc ){}; +}; + +/* Normal repulsion */ + +struct NormalRepulsionModel : public ContactModel +{ + using ContactModel::delta; + using ContactModel::Rc; + + double c; + double K; + + NormalRepulsionModel(){}; + NormalRepulsionModel( const double delta, const double Rc, const double _K ) + : ContactModel( delta, Rc ) + , K( _K ) + { + set_param( delta, Rc, K ); + } + + void set_param( const double _delta, const double _Rc, const double _K ) + { + delta = _delta; + Rc = _Rc; + K = _K; + // This could inherit from PMB (same c) + c = 18.0 * K / ( pi * delta * delta * delta * delta ); + } +}; + +/****************************************************************************** + Normal repulsion computation +******************************************************************************/ +template +class Force + : public Force +{ + public: + using base_type = Force; + + template + Force( const bool half_neigh, const ParticleType& particles, + const NormalRepulsionModel model ) + : base_type( half_neigh, model.Rc, particles.sliceCurrentPosition(), + particles.n_local, particles.ghost_mesh_lo, + particles.ghost_mesh_hi ) + , _model( model ) + { + for ( int d = 0; d < particles.dim; d++ ) + { + mesh_min[d] = particles.ghost_mesh_lo[d]; + mesh_max[d] = particles.ghost_mesh_hi[d]; + } + } + + template + void computeForceFull( 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.sliceVolume(); + const auto x = particles.sliceReferencePosition(); + const auto u = particles.sliceDisplacement(); + const auto y = particles.sliceCurrentPosition(); + + _neigh_list.build( y, 0, n_local, Rc, 1.0, mesh_min, mesh_max ); + + auto contact_full = KOKKOS_LAMBDA( const int i, const int j ) + { + double fcx_i = 0.0; + double fcy_i = 0.0; + double fcz_i = 0.0; + + double xi, r, s; + double rx, ry, rz; + getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz ); + + // Contact "stretch" + 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; + fcy_i = coeff * ry / r; + fcz_i = coeff * rz / r; + + fc( i, 0 ) += fcx_i; + fc( i, 1 ) += fcy_i; + fc( i, 2 ) += fcz_i; + }; + + // FIXME: using default space for now. + using exec_space = typename MemorySpace::execution_space; + Kokkos::RangePolicy policy( 0, n_local ); + Cabana::neighbor_parallel_for( + policy, contact_full, _neigh_list, Cabana::FirstNeighborsTag(), + neigh_op_tag, "CabanaPD::Contact::compute_full" ); + } + + protected: + NormalRepulsionModel _model; + using base_type::_half_neigh; + using base_type::_neigh_list; + + double mesh_max[3]; + double mesh_min[3]; +}; + +template +void computeContact( const ForceType& force, ParticleType& particles, + const ParallelType& neigh_op_tag ) +{ + auto n_local = particles.n_local; + auto f = particles.sliceForce(); + auto f_a = particles.sliceForceAtomic(); + + // 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 + // compute_force_half( f_a, x, u, neigh_list, n_local, + // neigh_op_tag ); + + // Forces only atomic if using team threading + if constexpr ( std::is_same::value ) + force.computeContactFull( f_a, particles, n_local, neigh_op_tag ); + else + force.computeContactFull( f, particles, n_local, neigh_op_tag ); + Kokkos::fence(); +} + +} // namespace CabanaPD + +#endif diff --git a/src/CabanaPD_Force.hpp b/src/CabanaPD_Force.hpp index 6dab7f08..168a5353 100644 --- a/src/CabanaPD_Force.hpp +++ b/src/CabanaPD_Force.hpp @@ -156,7 +156,20 @@ class Force { } - // Primary constructor: use positions and construct neighbors. + // General constructor (necessary for contact, but could be used by any + // force routine). + template + Force( const bool half_neigh, const double delta, + const PositionType& positions, const std::size_t num_local, + const double mesh_min[3], const double mesh_max[3], + const double tol = 1e-14 ) + : _half_neigh( half_neigh ) + , _neigh_list( neighbor_list_type( positions, 0, num_local, delta + tol, + 1.0, mesh_min, mesh_max ) ) + { + } + + // Constructor which stores existing neighbors. template Force( const bool half_neigh, const NeighborListType& neighbors ) : _half_neigh( half_neigh ) diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index d7a99321..40eef2ab 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -73,6 +73,7 @@ #include #include +#include #include #include #include @@ -92,7 +93,7 @@ class SolverBase }; template + class ForceModel, class ContactModel = NoContact> class SolverElastic { public: @@ -111,6 +112,8 @@ class SolverElastic // Optional module types. using heat_transfer_type = HeatTransfer; + using contact_type = Force; + using contact_model_type = ContactModel; SolverElastic( input_type _inputs, std::shared_ptr _particles, @@ -118,6 +121,27 @@ class SolverElastic : inputs( _inputs ) , particles( _particles ) , _init_time( 0.0 ) + { + setup( force_model ); + } + + SolverElastic( input_type _inputs, + std::shared_ptr _particles, + force_model_type force_model, + contact_model_type contact_model ) + : inputs( _inputs ) + , particles( _particles ) + , _init_time( 0.0 ) + { + setup( force_model ); + + _neighbor_timer.start(); + contact = std::make_shared( inputs["half_neigh"], + *particles, contact_model ); + _neighbor_timer.stop(); + } + + void setup( force_model_type force_model ) { num_steps = inputs["num_steps"]; output_frequency = inputs["output_frequency"]; @@ -135,7 +159,17 @@ class SolverElastic typename force_model_type::thermal_type>::value ) force_model.update( particles->sliceTemperature() ); + // Create heat transfer if needed. + if constexpr ( is_heat_transfer< + typename force_model_type::thermal_type>::value ) + { + thermal_subcycle_steps = inputs["thermal_subcycle_steps"]; + heat_transfer = std::make_shared( + inputs["half_neigh"], force->_neigh_list, force_model ); + } + _neighbor_timer.start(); + // This will either be PD or DEM forces. force = std::make_shared( inputs["half_neigh"], *particles, force_model ); _neighbor_timer.stop(); @@ -255,6 +289,9 @@ class SolverElastic // Compute internal forces. updateForce(); + if constexpr ( is_contact::value ) + computeContact( *contact, *particles, neigh_iter_tag{} ); + // Add force boundary condition. if ( boundary_condition.forceUpdate() ) boundary_condition.apply( exec_space(), *particles, step * dt ); @@ -287,6 +324,9 @@ class SolverElastic // Compute internal forces. updateForce(); + if constexpr ( is_contact::value ) + computeContact( *contact, *particles, neigh_iter_tag{} ); + if constexpr ( is_temperature_dependent< typename force_model_type::thermal_type>::value ) comm->gatherTemperature(); @@ -422,6 +462,8 @@ class SolverElastic std::shared_ptr force; // Optional modules. std::shared_ptr heat_transfer; + std::shared_ptr contact; + contact_model_type contact_model; // Output files. std::string output_file; @@ -437,13 +479,14 @@ class SolverElastic }; template + class ForceModel, class ContactModel = NoContact> class SolverFracture - : public SolverElastic + : public SolverElastic { public: - using base_type = - SolverElastic; + using base_type = SolverElastic; using exec_space = typename base_type::exec_space; using memory_space = typename base_type::memory_space; @@ -455,23 +498,21 @@ class SolverFracture using neigh_iter_tag = Cabana::SerialOpTag; using input_type = typename base_type::input_type; - template + using contact_model_type = ContactModel; + SolverFracture( input_type _inputs, std::shared_ptr _particles, - force_model_type force_model, PrenotchType prenotch ) + force_model_type force_model ) : base_type( _inputs, _particles, force_model ) { init_mu(); - - // Create prenotch. - prenotch.create( exec_space{}, mu, *particles, force->_neigh_list ); - _init_time += prenotch.time(); } SolverFracture( input_type _inputs, std::shared_ptr _particles, - force_model_type force_model ) - : base_type( _inputs, _particles, force_model ) + force_model_type force_model, + contact_model_type contact_model ) + : base_type( _inputs, _particles, force_model, contact_model ) { init_mu(); } @@ -488,6 +529,14 @@ class SolverFracture _init_timer.stop(); } + template + void init_prenotch( Prenotch prenotch ) + { + // Create prenotch. + prenotch.create( exec_space{}, mu, *particles, force->_neigh_list ); + _init_time += prenotch.time(); + } + void init( const bool initial_output = true ) { // Compute initial internal forces and energy. @@ -498,6 +547,14 @@ class SolverFracture particles->output( 0, 0.0, output_reference ); } + template + void init( Prenotch prenotch, + const bool initial_output = true ) + { + init_prenotch( prenotch ); + init( initial_output ); + } + template void init( BoundaryType boundary_condition, const bool initial_output = true ) @@ -522,6 +579,14 @@ class SolverFracture particles->output( 0, 0.0, output_reference ); } + template + void init( BoundaryType boundary_condition, Prenotch prenotch, + const bool initial_output = true ) + { + init_prenotch( prenotch ); + init( boundary_condition, initial_output ); + } + template void run( BoundaryType boundary_condition ) { @@ -549,6 +614,9 @@ class SolverFracture // Compute internal forces. updateForce(); + if constexpr ( is_contact::value ) + computeContact( *contact, *particles, neigh_iter_tag{} ); + // Add force boundary condition. if ( boundary_condition.forceUpdate() ) boundary_condition.apply( exec_space{}, *particles, step * dt ); @@ -585,6 +653,9 @@ class SolverFracture // Compute internal forces. updateForce(); + if constexpr ( is_contact::value ) + computeContact( *contact, *particles, neigh_iter_tag{} ); + // Integrate - velocity Verlet second half. integrator->finalHalfStep( *particles ); @@ -637,6 +708,7 @@ class SolverFracture protected: using base_type::comm; + using base_type::contact; using base_type::force; using base_type::inputs; using base_type::integrator; @@ -651,37 +723,29 @@ class SolverFracture using base_type::print; }; +// =============================================================== + template + class ForceModel, class... ContactModelType> auto createSolverElastic( InputsType inputs, std::shared_ptr particles, - ForceModel model ) -{ - return std::make_shared< - SolverElastic>( - inputs, particles, model ); -} - -template -auto createSolverFracture( InputsType inputs, - std::shared_ptr particles, - ForceModel model ) + ForceModel model, ContactModelType... contact_model ) { - return std::make_shared< - SolverFracture>( - inputs, particles, model ); + return std::make_shared>( + inputs, particles, model, contact_model... ); } template + class ForceModel, class... ContactModelType> auto createSolverFracture( InputsType inputs, std::shared_ptr particles, - ForceModel model, PrenotchType prenotch ) + ForceModel model, ContactModelType... contact_model ) { return std::make_shared< - SolverFracture>( - inputs, particles, model, prenotch ); + SolverFracture>( inputs, particles, model, + contact_model... ); } } // namespace CabanaPD diff --git a/src/CabanaPD_Types.hpp b/src/CabanaPD_Types.hpp index 58894e93..eb12a430 100644 --- a/src/CabanaPD_Types.hpp +++ b/src/CabanaPD_Types.hpp @@ -14,7 +14,7 @@ namespace CabanaPD { -// Mechanics types. +// Fracture tags. struct Elastic { }; @@ -22,7 +22,24 @@ struct Fracture { }; -// Thermal types. +// Contact and DEM (contact without PD) tags. + +struct Contact +{ +}; +struct NoContact +{ +}; +template +struct is_contact : public std::false_type +{ +}; +template <> +struct is_contact : public std::true_type +{ +}; + +// Thermal tags. struct TemperatureIndependent { using base_type = TemperatureIndependent; @@ -58,7 +75,7 @@ struct is_heat_transfer : public std::true_type { }; -// Model types. +// Force model tags. struct PMB { }; From f4acdbc080eda00e298aabea915dbc6b56b625fd Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Tue, 19 Nov 2024 13:32:39 -0500 Subject: [PATCH 02/23] Rearrange contact files --- src/CabanaPD.hpp | 3 +- src/force/CabanaPD_ForceModels_Contact.hpp | 68 +++++++++++++++++++ .../CabanaPD_Force_Contact.hpp} | 46 +------------ 3 files changed, 71 insertions(+), 46 deletions(-) create mode 100644 src/force/CabanaPD_ForceModels_Contact.hpp rename src/{CabanaPD_Contact.hpp => force/CabanaPD_Force_Contact.hpp} (81%) diff --git a/src/CabanaPD.hpp b/src/CabanaPD.hpp index 380ab056..4acf096d 100644 --- a/src/CabanaPD.hpp +++ b/src/CabanaPD.hpp @@ -16,7 +16,6 @@ #include #include #include -#include #include #include #include @@ -30,8 +29,10 @@ #include #include +#include #include #include +#include #include #include diff --git a/src/force/CabanaPD_ForceModels_Contact.hpp b/src/force/CabanaPD_ForceModels_Contact.hpp new file mode 100644 index 00000000..98176a53 --- /dev/null +++ b/src/force/CabanaPD_ForceModels_Contact.hpp @@ -0,0 +1,68 @@ +/**************************************************************************** + * 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 * + ****************************************************************************/ + +#ifndef CONTACTMODELS_H +#define CONTACTMODELS_H + +#include + +#include +#include +#include + +namespace CabanaPD +{ +/****************************************************************************** + Contact model +******************************************************************************/ +struct ContactModel +{ + double delta; + double Rc; + + ContactModel(){}; + // PD horizon + // Contact radius + ContactModel( const double _delta, const double _Rc ) + : delta( _delta ) + , Rc( _Rc ){}; +}; + +/* Normal repulsion */ + +struct NormalRepulsionModel : public ContactModel +{ + using ContactModel::delta; + using ContactModel::Rc; + + double c; + double K; + + NormalRepulsionModel(){}; + NormalRepulsionModel( const double delta, const double Rc, const double _K ) + : ContactModel( delta, Rc ) + , K( _K ) + { + set_param( delta, Rc, K ); + } + + void set_param( const double _delta, const double _Rc, const double _K ) + { + delta = _delta; + Rc = _Rc; + K = _K; + // This could inherit from PMB (same c) + c = 18.0 * K / ( pi * delta * delta * delta * delta ); + } +}; +} // namespace CabanaPD + +#endif diff --git a/src/CabanaPD_Contact.hpp b/src/force/CabanaPD_Force_Contact.hpp similarity index 81% rename from src/CabanaPD_Contact.hpp rename to src/force/CabanaPD_Force_Contact.hpp index 40cd81f8..8c26c108 100644 --- a/src/CabanaPD_Contact.hpp +++ b/src/force/CabanaPD_Force_Contact.hpp @@ -21,51 +21,7 @@ namespace CabanaPD { /****************************************************************************** - Contact model -******************************************************************************/ -struct ContactModel -{ - double delta; - double Rc; - - ContactModel(){}; - // PD horizon - // Contact radius - ContactModel( const double _delta, const double _Rc ) - : delta( _delta ) - , Rc( _Rc ){}; -}; - -/* Normal repulsion */ - -struct NormalRepulsionModel : public ContactModel -{ - using ContactModel::delta; - using ContactModel::Rc; - - double c; - double K; - - NormalRepulsionModel(){}; - NormalRepulsionModel( const double delta, const double Rc, const double _K ) - : ContactModel( delta, Rc ) - , K( _K ) - { - set_param( delta, Rc, K ); - } - - void set_param( const double _delta, const double _Rc, const double _K ) - { - delta = _delta; - Rc = _Rc; - K = _K; - // This could inherit from PMB (same c) - c = 18.0 * K / ( pi * delta * delta * delta * delta ); - } -}; - -/****************************************************************************** - Normal repulsion computation + Normal repulsion forces ******************************************************************************/ template class Force From 9146ebc449cfa07e020335d64919c38b5e895cff Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Tue, 19 Nov 2024 21:27:24 -0500 Subject: [PATCH 03/23] fixup: unify contact and force --- src/CabanaPD_Force.hpp | 10 ++++--- src/CabanaPD_Solver.hpp | 9 +++--- src/force/CabanaPD_Force_Contact.hpp | 41 ++++++++-------------------- 3 files changed, 22 insertions(+), 38 deletions(-) diff --git a/src/CabanaPD_Force.hpp b/src/CabanaPD_Force.hpp index 168a5353..8299c78b 100644 --- a/src/CabanaPD_Force.hpp +++ b/src/CabanaPD_Force.hpp @@ -237,7 +237,7 @@ class Force ******************************************************************************/ template void computeForce( ForceType& force, ParticleType& particles, - const ParallelType& neigh_op_tag ) + const ParallelType& neigh_op_tag, const bool reset = true ) { auto n_local = particles.n_local; auto x = particles.sliceReferencePosition(); @@ -246,7 +246,8 @@ void computeForce( ForceType& force, ParticleType& particles, auto f_a = particles.sliceForceAtomic(); // Reset force. - Cabana::deep_copy( f, 0.0 ); + if ( reset ) + Cabana::deep_copy( f, 0.0 ); // if ( half_neigh ) // Forces must be atomic for half list @@ -291,7 +292,7 @@ double computeEnergy( ForceType& force, ParticleType& particles, template void computeForce( ForceType& force, ParticleType& particles, NeighborView& mu, - const ParallelType& neigh_op_tag ) + const ParallelType& neigh_op_tag, const bool reset = true ) { auto n_local = particles.n_local; auto x = particles.sliceReferencePosition(); @@ -300,7 +301,8 @@ void computeForce( ForceType& force, ParticleType& particles, NeighborView& mu, auto f_a = particles.sliceForceAtomic(); // Reset force. - Cabana::deep_copy( f, 0.0 ); + if ( reset ) + Cabana::deep_copy( f, 0.0 ); // if ( half_neigh ) // Forces must be atomic for half list diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 40eef2ab..a1db89a2 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -73,7 +73,6 @@ #include #include -#include #include #include #include @@ -290,7 +289,7 @@ class SolverElastic updateForce(); if constexpr ( is_contact::value ) - computeContact( *contact, *particles, neigh_iter_tag{} ); + computeForce( *contact, *particles, neigh_iter_tag{}, false ); // Add force boundary condition. if ( boundary_condition.forceUpdate() ) @@ -325,7 +324,7 @@ class SolverElastic updateForce(); if constexpr ( is_contact::value ) - computeContact( *contact, *particles, neigh_iter_tag{} ); + computeForce( *contact, *particles, neigh_iter_tag{}, false ); if constexpr ( is_temperature_dependent< typename force_model_type::thermal_type>::value ) @@ -615,7 +614,7 @@ class SolverFracture updateForce(); if constexpr ( is_contact::value ) - computeContact( *contact, *particles, neigh_iter_tag{} ); + computeForce( *contact, *particles, neigh_iter_tag{}, false ); // Add force boundary condition. if ( boundary_condition.forceUpdate() ) @@ -654,7 +653,7 @@ class SolverFracture updateForce(); if constexpr ( is_contact::value ) - computeContact( *contact, *particles, neigh_iter_tag{} ); + computeForce( *contact, *particles, neigh_iter_tag{}, false ); // Integrate - velocity Verlet second half. integrator->finalHalfStep( *particles ); diff --git a/src/force/CabanaPD_Force_Contact.hpp b/src/force/CabanaPD_Force_Contact.hpp index 8c26c108..628c4c09 100644 --- a/src/force/CabanaPD_Force_Contact.hpp +++ b/src/force/CabanaPD_Force_Contact.hpp @@ -45,19 +45,21 @@ class Force } } - template + template void computeForceFull( ForceType& fc, const ParticleType& particles, - const int n_local, ParallelType& neigh_op_tag ) const + const PosType& x, const PosType& u, + const int n_local, ParallelType& neigh_op_tag ) { auto delta = _model.delta; auto Rc = _model.Rc; auto c = _model.c; const auto vol = particles.sliceVolume(); - const auto x = particles.sliceReferencePosition(); - const auto u = particles.sliceDisplacement(); const auto y = particles.sliceCurrentPosition(); + _neigh_timer.start(); _neigh_list.build( y, 0, n_local, Rc, 1.0, mesh_min, mesh_max ); + _neigh_timer.stop(); auto contact_full = KOKKOS_LAMBDA( const int i, const int j ) { @@ -83,48 +85,29 @@ class Force fc( i, 2 ) += fcz_i; }; + _timer.start(); + // FIXME: using default space for now. using exec_space = typename MemorySpace::execution_space; Kokkos::RangePolicy policy( 0, n_local ); Cabana::neighbor_parallel_for( policy, contact_full, _neigh_list, Cabana::FirstNeighborsTag(), neigh_op_tag, "CabanaPD::Contact::compute_full" ); + + _timer.stop(); } protected: NormalRepulsionModel _model; using base_type::_half_neigh; using base_type::_neigh_list; + using base_type::_timer; + Timer _neigh_timer; double mesh_max[3]; double mesh_min[3]; }; -template -void computeContact( const ForceType& force, ParticleType& particles, - const ParallelType& neigh_op_tag ) -{ - auto n_local = particles.n_local; - auto f = particles.sliceForce(); - auto f_a = particles.sliceForceAtomic(); - - // 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 - // compute_force_half( f_a, x, u, neigh_list, n_local, - // neigh_op_tag ); - - // Forces only atomic if using team threading - if constexpr ( std::is_same::value ) - force.computeContactFull( f_a, particles, n_local, neigh_op_tag ); - else - force.computeContactFull( f, particles, n_local, neigh_op_tag ); - Kokkos::fence(); -} - } // namespace CabanaPD #endif From 4d039bec8368ab7a8f50e18e8ebba74e0c682f3a Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 20 Nov 2024 09:39:51 -0500 Subject: [PATCH 04/23] Fix terrible mistake in impact inputs --- examples/mechanics/inputs/disk_impact.json | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/mechanics/inputs/disk_impact.json b/examples/mechanics/inputs/disk_impact.json index 6a273f0b..b02a77c1 100644 --- a/examples/mechanics/inputs/disk_impact.json +++ b/examples/mechanics/inputs/disk_impact.json @@ -1,6 +1,6 @@ { "num_cells" : {"value": [149, 149, 5]}, - "low_corner" : {"value": [-0.037, -0.037, 2.5e-3], "unit": "m"}, + "low_corner" : {"value": [-0.037, -0.037, -2.5e-3], "unit": "m"}, "high_corner" : {"value": [0.037, 0.037, 0.0], "unit": "m"}, "density" : {"value": 2200.0, "unit": "kg/m^3"}, "bulk_modulus" : {"value": 14.9e+9, "unit": "Pa"}, From 340987b5824178fbcd8eeed9a5a506a5e9c2e368 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 20 Nov 2024 09:40:10 -0500 Subject: [PATCH 05/23] Support system size and corners inputs --- src/CabanaPD_Input.hpp | 67 +++++++++++++++++++++++++++++++----------- 1 file changed, 50 insertions(+), 17 deletions(-) diff --git a/src/CabanaPD_Input.hpp b/src/CabanaPD_Input.hpp index b615b15e..596d42ee 100644 --- a/src/CabanaPD_Input.hpp +++ b/src/CabanaPD_Input.hpp @@ -32,22 +32,7 @@ class Inputs inputs = parse( filename ); // Add additional derived inputs to json. System size. - auto size = inputs["system_size"]["value"]; - std::string size_unit = inputs["system_size"]["unit"]; - for ( std::size_t d = 0; d < size.size(); d++ ) - { - double s = size[d]; - double low = -0.5 * s; - double high = 0.5 * s; - inputs["low_corner"]["value"][d] = low; - inputs["high_corner"]["value"][d] = high; - - double nc = inputs["num_cells"]["value"][d]; - inputs["dx"]["value"][d] = ( high - low ) / nc; - } - inputs["low_corner"]["unit"] = size_unit; - inputs["high_corner"]["unit"] = size_unit; - inputs["dx"]["unit"] = size_unit; + setupSize(); // Number of steps. double tf = inputs["final_time"]["value"]; @@ -104,7 +89,55 @@ class Inputs // Not yet a user option. inputs["half_neigh"]["value"] = false; } - ~Inputs() {} + + void setupSize() + { + if ( inputs.contains( "system_size" ) ) + { + auto size = inputs["system_size"]["value"]; + std::string size_unit = inputs["system_size"]["unit"]; + for ( std::size_t d = 0; d < size.size(); d++ ) + { + double s = size[d]; + double low = -0.5 * s; + double high = 0.5 * s; + inputs["low_corner"]["value"][d] = low; + inputs["high_corner"]["value"][d] = high; + } + inputs["low_corner"]["unit"] = size_unit; + inputs["high_corner"]["unit"] = size_unit; + } + else if ( inputs.contains( "low_corner" ) && + ( inputs.contains( "high_corner" ) ) ) + { + auto low_corner = inputs["low_corner"]["value"]; + auto high_corner = inputs["high_corner"]["value"]; + std::string size_unit = inputs["low_corner"]["unit"]; + for ( std::size_t d = 0; d < low_corner.size(); d++ ) + { + double low = low_corner[d]; + double high = high_corner[d]; + inputs["system_size"]["value"][d] = high - low; + } + inputs["system_size"]["unit"] = size_unit; + } + else + { + throw std::runtime_error( "Must input either system_size or " + "both low_corner and high_corner." ); + } + + auto size = inputs["system_size"]["value"]; + for ( std::size_t d = 0; d < size.size(); d++ ) + { + double low = inputs["low_corner"]["value"][d]; + double high = inputs["low_corner"]["value"][d]; + double nc = inputs["num_cells"]["value"][d]; + inputs["dx"]["value"][d] = ( high - low ) / nc; + } + std::string size_unit = inputs["system_size"]["unit"]; + inputs["dx"]["unit"] = size_unit; + } void computeCriticalTimeStep() { From d6e4ad97c508ef2af991431e2a0fa5b6bd8792c0 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 20 Nov 2024 09:40:33 -0500 Subject: [PATCH 06/23] cleanup impact problem --- examples/mechanics/disk_impact.cpp | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/examples/mechanics/disk_impact.cpp b/examples/mechanics/disk_impact.cpp index 0877b4a2..42fee4c3 100644 --- a/examples/mechanics/disk_impact.cpp +++ b/examples/mechanics/disk_impact.cpp @@ -28,15 +28,12 @@ void diskImpactExample( const std::string filename ) auto particles = CabanaPD::createParticles( exec_space(), inputs ); - std::array low_corner = inputs["low_corner"]; - std::array high_corner = inputs["high_corner"]; - double global_mesh_ext = high_corner[0] - low_corner[0]; - auto x = particles->sliceReferencePosition(); + double system_size = inputs["system_size"][0]; auto create_func = KOKKOS_LAMBDA( const int, const double px[3] ) { - auto width = global_mesh_ext / 2.0; + auto radius = system_size / 2.0 - 1e-10; auto r2 = px[0] * px[0] + px[1] * px[1]; - if ( r2 > width * width ) + if ( r2 > radius * radius ) return false; return true; }; @@ -44,9 +41,7 @@ void diskImpactExample( const std::string filename ) // Define particle initialization. auto v = particles->sliceVelocity(); - auto f = particles->sliceForce(); auto rho = particles->sliceDensity(); - double rho0 = inputs["density"]; auto init_functor = KOKKOS_LAMBDA( const int pid ) { @@ -60,7 +55,6 @@ void diskImpactExample( const std::string filename ) double r_c = dx * 2.0; CabanaPD::NormalRepulsionModel contact_model( delta, r_c, K ); - // FIXME: use createSolver to switch backend at runtime. auto cabana_pd = CabanaPD::createSolverFracture( inputs, particles, force_model, contact_model ); @@ -68,6 +62,8 @@ void diskImpactExample( const std::string filename ) double impact_v = inputs["impactor_velocity"]; double impact_x = 0.0; double impact_y = 0.0; + auto x = particles->sliceReferencePosition(); + auto f = particles->sliceForce(); auto body_func = KOKKOS_LAMBDA( const int p, const double t ) { auto z = t * impact_v + impact_r; From 9967b8ab07e7e3409b4f9257754bc22e84e7ebb1 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 20 Nov 2024 09:43:17 -0500 Subject: [PATCH 07/23] Add contact to fragmentation --- examples/mechanics/fragmenting_cylinder.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/examples/mechanics/fragmenting_cylinder.cpp b/examples/mechanics/fragmenting_cylinder.cpp index 0c84a553..083e954c 100644 --- a/examples/mechanics/fragmenting_cylinder.cpp +++ b/examples/mechanics/fragmenting_cylinder.cpp @@ -136,11 +136,14 @@ void fragmentingCylinderExample( const std::string filename ) }; particles->updateParticles( exec_space{}, init_functor ); + double r_c = dx[0] * 2.0; + CabanaPD::NormalRepulsionModel contact_model( delta, r_c, K ); + // ==================================================== // Simulation run // ==================================================== auto cabana_pd = CabanaPD::createSolverFracture( - inputs, particles, force_model ); + inputs, particles, force_model, contact_model ); cabana_pd->init(); cabana_pd->run(); } From 8fd51669b7e0eb2d4927c8c408b350ca9a6ab4b2 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 20 Nov 2024 16:38:55 -0500 Subject: [PATCH 08/23] Support dx or num_cells input --- src/CabanaPD_Input.hpp | 34 ++++++++++++++++++++++++++-------- 1 file changed, 26 insertions(+), 8 deletions(-) diff --git a/src/CabanaPD_Input.hpp b/src/CabanaPD_Input.hpp index 596d42ee..29a0f898 100644 --- a/src/CabanaPD_Input.hpp +++ b/src/CabanaPD_Input.hpp @@ -127,16 +127,34 @@ class Inputs "both low_corner and high_corner." ); } - auto size = inputs["system_size"]["value"]; - for ( std::size_t d = 0; d < size.size(); d++ ) + if ( inputs.contains( "dx" ) ) { - double low = inputs["low_corner"]["value"][d]; - double high = inputs["low_corner"]["value"][d]; - double nc = inputs["num_cells"]["value"][d]; - inputs["dx"]["value"][d] = ( high - low ) / nc; + auto size = inputs["system_size"]["value"]; + for ( std::size_t d = 0; d < size.size(); d++ ) + { + double system_size = inputs["system_size"]["value"][d]; + double dx = inputs["dx"]["value"][d]; + inputs["num_cells"]["value"][d] = + static_cast( system_size / dx ); + } + } + else if ( inputs.contains( "num_cells" ) ) + { + auto size = inputs["system_size"]["value"]; + for ( std::size_t d = 0; d < size.size(); d++ ) + { + double low = inputs["low_corner"]["value"][d]; + double high = inputs["low_corner"]["value"][d]; + double nc = inputs["num_cells"]["value"][d]; + inputs["dx"]["value"][d] = ( high - low ) / nc; + } + std::string size_unit = inputs["system_size"]["unit"]; + inputs["dx"]["unit"] = size_unit; + } + else + { + throw std::runtime_error( "Must input either num_cells or dx." ); } - std::string size_unit = inputs["system_size"]["unit"]; - inputs["dx"]["unit"] = size_unit; } void computeCriticalTimeStep() From ba1e048900efdeb656dba3ef31a2745a058c4692 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 20 Nov 2024 16:40:44 -0500 Subject: [PATCH 09/23] Fix contact static checking; associated issues in using contact at all --- src/CabanaPD_Types.hpp | 7 ------- src/force/CabanaPD_ForceModels_Contact.hpp | 11 +++++++++++ src/force/CabanaPD_Force_Contact.hpp | 15 ++++++++++++--- 3 files changed, 23 insertions(+), 10 deletions(-) diff --git a/src/CabanaPD_Types.hpp b/src/CabanaPD_Types.hpp index eb12a430..4ca10d35 100644 --- a/src/CabanaPD_Types.hpp +++ b/src/CabanaPD_Types.hpp @@ -24,9 +24,6 @@ struct Fracture // Contact and DEM (contact without PD) tags. -struct Contact -{ -}; struct NoContact { }; @@ -34,10 +31,6 @@ template struct is_contact : public std::false_type { }; -template <> -struct is_contact : public std::true_type -{ -}; // Thermal tags. struct TemperatureIndependent diff --git a/src/force/CabanaPD_ForceModels_Contact.hpp b/src/force/CabanaPD_ForceModels_Contact.hpp index 98176a53..b588e52d 100644 --- a/src/force/CabanaPD_ForceModels_Contact.hpp +++ b/src/force/CabanaPD_ForceModels_Contact.hpp @@ -40,6 +40,11 @@ struct ContactModel struct NormalRepulsionModel : public ContactModel { + // FIXME: This is for use as the primary force model. + using base_model = PMB; + using fracture_type = Elastic; + using thermal_type = TemperatureIndependent; + using ContactModel::delta; using ContactModel::Rc; @@ -63,6 +68,12 @@ struct NormalRepulsionModel : public ContactModel c = 18.0 * K / ( pi * delta * delta * delta * delta ); } }; + +template <> +struct is_contact : public std::true_type +{ +}; + } // namespace CabanaPD #endif diff --git a/src/force/CabanaPD_Force_Contact.hpp b/src/force/CabanaPD_Force_Contact.hpp index 628c4c09..b04956b3 100644 --- a/src/force/CabanaPD_Force_Contact.hpp +++ b/src/force/CabanaPD_Force_Contact.hpp @@ -29,6 +29,7 @@ class Force { public: using base_type = Force; + using neighbor_list_type = typename base_type::neighbor_list_type; template Force( const bool half_neigh, const ParticleType& particles, @@ -47,9 +48,9 @@ class Force template - void computeForceFull( ForceType& fc, const ParticleType& particles, - const PosType& x, const PosType& u, - const int n_local, ParallelType& neigh_op_tag ) + void computeForceFull( ForceType& fc, const PosType& x, const PosType& u, + const ParticleType& particles, const int n_local, + ParallelType& neigh_op_tag ) { auto delta = _model.delta; auto Rc = _model.Rc; @@ -97,6 +98,14 @@ class Force _timer.stop(); } + // FIXME: implement energy + template + double computeEnergyFull( WType&, const PosType&, const PosType&, + ParticleType&, const int, ParallelType& ) + { + } + protected: NormalRepulsionModel _model; using base_type::_half_neigh; From 4a9e0b984877fa3ddc095f7bb8a8ef6326f09ec2 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 20 Nov 2024 16:41:15 -0500 Subject: [PATCH 10/23] fixup: inputs for problems with contact --- examples/mechanics/disk_impact.cpp | 8 ++++---- examples/mechanics/fragmenting_cylinder.cpp | 3 ++- examples/mechanics/inputs/disk_impact.json | 9 +++++---- examples/mechanics/inputs/fragmenting_cylinder.json | 1 + 4 files changed, 12 insertions(+), 9 deletions(-) diff --git a/examples/mechanics/disk_impact.cpp b/examples/mechanics/disk_impact.cpp index 42fee4c3..bef00871 100644 --- a/examples/mechanics/disk_impact.cpp +++ b/examples/mechanics/disk_impact.cpp @@ -33,7 +33,7 @@ void diskImpactExample( const std::string filename ) { auto radius = system_size / 2.0 - 1e-10; auto r2 = px[0] * px[0] + px[1] * px[1]; - if ( r2 > radius * radius ) + if ( r2 > radius * radius || px[2] < 0.0 ) return false; return true; }; @@ -51,8 +51,8 @@ void diskImpactExample( const std::string filename ) }; particles->updateParticles( exec_space{}, init_functor ); - double dx = particles->dx[0]; - double r_c = dx * 2.0; + double r_c = inputs["contact_horizon_factor"]; + r_c *= particles->dx[0]; CabanaPD::NormalRepulsionModel contact_model( delta, r_c, K ); auto cabana_pd = CabanaPD::createSolverFracture( @@ -72,7 +72,7 @@ void diskImpactExample( const std::string filename ) ( x( p, 2 ) - z ) * ( x( p, 2 ) - z ) ); if ( r < impact_r ) { - double fmag = -1.0e17 * ( r - impact_r ) * ( r - impact_r ); + double fmag = 1.0e17 * ( r - impact_r ) * ( r - impact_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; diff --git a/examples/mechanics/fragmenting_cylinder.cpp b/examples/mechanics/fragmenting_cylinder.cpp index 083e954c..ea190e9a 100644 --- a/examples/mechanics/fragmenting_cylinder.cpp +++ b/examples/mechanics/fragmenting_cylinder.cpp @@ -136,7 +136,8 @@ void fragmentingCylinderExample( const std::string filename ) }; particles->updateParticles( exec_space{}, init_functor ); - double r_c = dx[0] * 2.0; + double r_c = inputs["contact_horizon_factor"]; + r_c *= dx[0]; CabanaPD::NormalRepulsionModel contact_model( delta, r_c, K ); // ==================================================== diff --git a/examples/mechanics/inputs/disk_impact.json b/examples/mechanics/inputs/disk_impact.json index b02a77c1..dc2305e7 100644 --- a/examples/mechanics/inputs/disk_impact.json +++ b/examples/mechanics/inputs/disk_impact.json @@ -1,15 +1,16 @@ { - "num_cells" : {"value": [149, 149, 5]}, - "low_corner" : {"value": [-0.037, -0.037, -2.5e-3], "unit": "m"}, - "high_corner" : {"value": [0.037, 0.037, 0.0], "unit": "m"}, + "dx" : {"value": [0.0005, 0.0005, 0.0005], "note": "This is used to enable a larger box in Z."}, + "low_corner" : {"value": [-0.03725, -0.03725, -10.0e-2], "unit": "m"}, + "high_corner" : {"value": [0.03725, 0.03725, 2.5e-3], "unit": "m"}, "density" : {"value": 2200.0, "unit": "kg/m^3"}, "bulk_modulus" : {"value": 14.9e+9, "unit": "Pa"}, "fracture_energy" : {"value": 10.0575, "unit": "J/m^2"}, "horizon" : {"value": 0.0015, "unit": "m"}, + "contact_horizon_factor" : {"value": 0.9}, "final_time" : {"value": 2.0e-4, "unit": "s"}, "timestep" : {"value": 1.0e-7, "unit": "s"}, "timestep_safety_factor" : {"value": 0.85}, - "output_frequency" : {"value": 50}, + "output_frequency" : {"value": 25}, "output_reference" : {"value": false}, "impactor_radius" : {"value": 5e-3, "unit": "m"}, "impactor_velocity" : {"value": -100.0, "unit": "m/s"} diff --git a/examples/mechanics/inputs/fragmenting_cylinder.json b/examples/mechanics/inputs/fragmenting_cylinder.json index a6b7041b..e2c945b3 100644 --- a/examples/mechanics/inputs/fragmenting_cylinder.json +++ b/examples/mechanics/inputs/fragmenting_cylinder.json @@ -7,6 +7,7 @@ "shear_modulus" : {"value": 78e+9, "unit": "Pa"}, "critical_stretch" : {"value": 0.02}, "horizon" : {"value": 0.00417462, "unit": "m"}, + "contact_horizon_factor" : {"value": 0.9}, "cylinder_outer_radius" : {"value": 0.025, "unit": "m"}, "cylinder_inner_radius" : {"value": 0.02, "unit": "m"}, "max_radial_velocity" : {"value": 200, "unit": "m/s"}, From e62df9482cbde20376427809569fb88835dce7d0 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 20 Nov 2024 17:03:25 -0500 Subject: [PATCH 11/23] fixup: missing return --- src/force/CabanaPD_Force_Contact.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/force/CabanaPD_Force_Contact.hpp b/src/force/CabanaPD_Force_Contact.hpp index b04956b3..f3b389ef 100644 --- a/src/force/CabanaPD_Force_Contact.hpp +++ b/src/force/CabanaPD_Force_Contact.hpp @@ -104,6 +104,7 @@ class Force double computeEnergyFull( WType&, const PosType&, const PosType&, ParticleType&, const int, ParallelType& ) { + return 0.0; } protected: From e4eb7650d3b02da67f6a7865422bad9c32a9e5f8 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 25 Nov 2024 08:40:35 -0500 Subject: [PATCH 12/23] WIP add contact comm --- src/CabanaPD_Comm.hpp | 92 +++++++++++++++++++++++++++++++------- src/CabanaPD_Particles.hpp | 8 ++-- src/CabanaPD_Solver.hpp | 25 +++++++++++ src/CabanaPD_Types.hpp | 3 ++ 4 files changed, 109 insertions(+), 19 deletions(-) diff --git a/src/CabanaPD_Comm.hpp b/src/CabanaPD_Comm.hpp index 6ab1508c..9e8d1654 100644 --- a/src/CabanaPD_Comm.hpp +++ b/src/CabanaPD_Comm.hpp @@ -195,6 +195,9 @@ struct HaloIds Kokkos::parallel_for( "CabanaPD::Comm::GhostSearch", policy, ghost_search ); Kokkos::fence(); + + // Rebuild if needed. + rebuild( positions ); } template @@ -243,6 +246,7 @@ class Comm int max_export; using memory_space = typename ParticleType::memory_space; + using local_grid_type = typename ParticleType::local_grid_type; using halo_type = Cabana::Halo; using gather_u_type = Cabana::Gather; @@ -262,12 +266,11 @@ class Comm auto halo_width = local_grid->haloCellWidth(); auto topology = Cabana::Grid::getTopology( *local_grid ); - // Determine which particles need to be ghosted to neighbors. + // Determine which particles should be ghosted, reallocating and + // recounting if needed. // FIXME: set halo width based on cutoff distance. - auto halo_ids = - createHaloIds( *local_grid, positions, halo_width, max_export ); - // Rebuild if needed. - halo_ids.rebuild( positions ); + HaloIds halo_ids( + *local_grid, positions, halo_width, max_export ); // Create the Cabana Halo. halo = std::make_shared( local_grid->globalGrid().comm(), @@ -288,17 +291,6 @@ class Comm } ~Comm() {} - // Determine which particles should be ghosted, reallocating and recounting - // if needed. - template - auto createHaloIds( const LocalGridType& local_grid, - const PositionSliceType& positions, - const int min_halo_width, const int max_export ) - { - return HaloIds( - local_grid, positions, min_halo_width, max_export ); - } - // We assume here that the particle count has not changed and no resize // is necessary. void gatherDisplacement() @@ -393,6 +385,74 @@ class Comm void gatherTemperature() { gather_temp->apply(); } }; +// Does not inherit because it does not use the same reference halo +// communication pattern. +template +class Comm +{ + public: + using memory_space = typename ParticleType::memory_space; + using local_grid_type = typename ParticleType::local_grid_type; + using halo_type = Cabana::Halo; + + using gather_current_type = + Cabana::Gather; + std::shared_ptr gather_current; + std::shared_ptr halo; + HaloIds halo_ids; + int max_export; + int halo_width; + + // Note this initial guess is small because this is often used for very + // short range interactions. + Comm( ParticleType& particles, int max_export_guess = 10 ) + : halo_ids( HaloIds( + *( particles.local_grid ), particles.sliceCurrentPosition(), + particles.local_grid->haloCellWidth(), max_export_guess ) ) + , max_export( max_export_guess ) + , halo_width( particles.local_grid->haloCellWidth() ) + { + // We use n_ghost here as the "local" halo count because these current + // frame ghosts are built on top of the existing, static, reference + // frame ghosts. + auto n_reference = particles.n_local + particles.n_ghost; + auto topology = Cabana::Grid::getTopology( *particles.local_grid ); + halo = std::make_shared( + particles.local_grid->globalGrid().comm(), n_reference, + halo_ids._ids, halo_ids._destinations, topology ); + std::cout << "halo " << halo->numLocal() << " " << halo->numGhost() + << std::endl; + particles.resize( halo->numLocal(), halo->numGhost() ); + + gather_current = + std::make_shared( *halo, particles._aosoa_y ); + } + + // This is a dynamic gather step where the steering vector needs to be + // recomputed. + void gatherCurrentPostiion( ParticleType& particles ) + { + // Get the current position. Note this is necessary to get the up to + // date current position. + auto y = particles.sliceCurrentPosition(); + // Determine which particles need to be ghosted to neighbors for the + // current positions. + halo_ids.build( y ); + + auto n_reference = particles.n_local + particles.n_ghost; + auto topology = Cabana::Grid::getTopology( *particles.local_grid ); + // FIXME: missing a build() interface + halo = std::make_shared( + particles.local_grid->globalGrid().comm(), n_reference, + halo_ids._ids, halo_ids._destinations, topology ); + particles.resize( halo->numLocal(), halo->numGhost() ); + + gather_current->reserve( *halo, particles._aosoa_y ); + + gather_current->apply(); + } +}; + } // namespace CabanaPD #endif diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 52cc36ed..b3db1847 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -137,9 +137,10 @@ class Particles // FIXME: this is for neighborlist construction. double ghost_mesh_lo[dim]; double ghost_mesh_hi[dim]; - std::shared_ptr< - Cabana::Grid::LocalGrid>> - local_grid; + + using local_grid_type = + Cabana::Grid::LocalGrid>; + std::shared_ptr local_grid; Kokkos::Array dx; int halo_width; @@ -476,6 +477,7 @@ class Particles friend class Comm; friend class Comm; + friend class Comm; protected: aosoa_u_type _aosoa_u; diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index a1db89a2..5f30dd38 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -113,6 +113,8 @@ class SolverElastic using heat_transfer_type = HeatTransfer; using contact_type = Force; using contact_model_type = ContactModel; + using contact_comm_type = + Comm; SolverElastic( input_type _inputs, std::shared_ptr _particles, @@ -150,9 +152,20 @@ class SolverElastic dt = inputs["timestep"]; integrator = std::make_shared( dt ); + std::cout << particles->n_local << " " << particles->n_ghost + << std::endl; // Add ghosts from other MPI ranks. comm = std::make_shared( *particles ); + std::cout << particles->n_local << " " << particles->n_ghost + << std::endl; + + if constexpr ( is_contact::value ) + contact_comm = std::make_shared( *particles ); + + std::cout << particles->n_local << " " << particles->n_ghost + << std::endl; + // Update temperature ghost size if needed. if constexpr ( is_temperature_dependent< typename force_model_type::thermal_type>::value ) @@ -462,6 +475,7 @@ class SolverElastic // Optional modules. std::shared_ptr heat_transfer; std::shared_ptr contact; + std::shared_ptr contact_comm; contact_model_type contact_model; // Output files. @@ -610,6 +624,11 @@ class SolverFracture // Update ghost particles. comm->gatherDisplacement(); + // Need to ghost current positions for DEM/contact. + if constexpr ( is_contact::value || + is_contact::value ) + contact_comm->gatherCurrentPostiion( *particles ); + // Compute internal forces. updateForce(); @@ -649,6 +668,11 @@ class SolverFracture // Update ghost particles. comm->gatherDisplacement(); + // Need to ghost current positions for DEM/contact. + if constexpr ( is_contact::value || + is_contact::value ) + contact_comm->gatherCurrentPostiion( *particles ); + // Compute internal forces. updateForce(); @@ -708,6 +732,7 @@ class SolverFracture protected: using base_type::comm; using base_type::contact; + using base_type::contact_comm; using base_type::force; using base_type::inputs; using base_type::integrator; diff --git a/src/CabanaPD_Types.hpp b/src/CabanaPD_Types.hpp index 4ca10d35..2e2c5c98 100644 --- a/src/CabanaPD_Types.hpp +++ b/src/CabanaPD_Types.hpp @@ -24,6 +24,9 @@ struct Fracture // Contact and DEM (contact without PD) tags. +struct Contact +{ +}; struct NoContact { }; From c4bc765806e9d45668778939bb7e4ba1e13b66bc Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 25 Nov 2024 08:40:52 -0500 Subject: [PATCH 13/23] fixup: fragmentation box inputs --- examples/mechanics/fragmenting_cylinder.cpp | 2 +- examples/mechanics/inputs/fragmenting_cylinder.json | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/mechanics/fragmenting_cylinder.cpp b/examples/mechanics/fragmenting_cylinder.cpp index ea190e9a..64c6dcae 100644 --- a/examples/mechanics/fragmenting_cylinder.cpp +++ b/examples/mechanics/fragmenting_cylinder.cpp @@ -85,7 +85,7 @@ void fragmentingCylinderExample( const std::string filename ) { double rsq = ( x[0] - x_center ) * ( x[0] - x_center ) + ( x[1] - y_center ) * ( x[1] - y_center ); - if ( rsq < Rin * Rin || rsq > Rout * Rout ) + if ( rsq < Rin * Rin || rsq > Rout * Rout || x[2] > 0.05 || x[2] < -0.05) return false; return true; }; diff --git a/examples/mechanics/inputs/fragmenting_cylinder.json b/examples/mechanics/inputs/fragmenting_cylinder.json index e2c945b3..b33b28f5 100644 --- a/examples/mechanics/inputs/fragmenting_cylinder.json +++ b/examples/mechanics/inputs/fragmenting_cylinder.json @@ -1,6 +1,7 @@ { "num_cells" : {"value": [51, 51, 101]}, - "system_size" : {"value": [0.05, 0.05, 0.1], "unit": "m"}, + "dx": {"value": [0.001, 0.001, 0.001]}, + "system_size" : {"value": [0.1, 0.1, 0.15], "unit": "m"}, "grid_perturbation_factor" : {"value": 0.1}, "density" : {"value": 7800, "unit": "kg/m^3"}, "bulk_modulus" : {"value": 130e+9, "unit": "Pa"}, From 8ebf293b1f09caac7d54afa2f4ded3e3572977f7 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 27 Nov 2024 10:49:00 -0600 Subject: [PATCH 14/23] fixup: local, ghost, and contact ghost clarifications --- src/CabanaPD_Comm.hpp | 18 +++++++++--------- src/CabanaPD_Particles.hpp | 35 ++++++++++++++++++++++------------- src/CabanaPD_Solver.hpp | 14 ++++++++------ 3 files changed, 39 insertions(+), 28 deletions(-) diff --git a/src/CabanaPD_Comm.hpp b/src/CabanaPD_Comm.hpp index 9e8d1654..22c2686c 100644 --- a/src/CabanaPD_Comm.hpp +++ b/src/CabanaPD_Comm.hpp @@ -412,17 +412,17 @@ class Comm , max_export( max_export_guess ) , halo_width( particles.local_grid->haloCellWidth() ) { - // We use n_ghost here as the "local" halo count because these current - // frame ghosts are built on top of the existing, static, reference - // frame ghosts. - auto n_reference = particles.n_local + particles.n_ghost; auto topology = Cabana::Grid::getTopology( *particles.local_grid ); halo = std::make_shared( - particles.local_grid->globalGrid().comm(), n_reference, + particles.local_grid->globalGrid().comm(), particles.n_reference, halo_ids._ids, halo_ids._destinations, topology ); std::cout << "halo " << halo->numLocal() << " " << halo->numGhost() << std::endl; - particles.resize( halo->numLocal(), halo->numGhost() ); + // We use n_ghost here as the "local" halo count because these current + // frame ghosts are built on top of the existing, static, reference + // frame ghosts. + particles.resize( particles.n_local, particles.n_ghost, + halo->numGhost() ); gather_current = std::make_shared( *halo, particles._aosoa_y ); @@ -439,13 +439,13 @@ class Comm // current positions. halo_ids.build( y ); - auto n_reference = particles.n_local + particles.n_ghost; auto topology = Cabana::Grid::getTopology( *particles.local_grid ); // FIXME: missing a build() interface halo = std::make_shared( - particles.local_grid->globalGrid().comm(), n_reference, + particles.local_grid->globalGrid().comm(), particles.n_reference, halo_ids._ids, halo_ids._destinations, topology ); - particles.resize( halo->numLocal(), halo->numGhost() ); + particles.resize( particles.n_local, particles.n_ghost, + halo->numGhost() ); gather_current->reserve( *halo, particles._aosoa_y ); diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index b3db1847..24844aca 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -97,6 +97,8 @@ class Particles unsigned long long int n_global = 0; std::size_t n_local = 0; std::size_t n_ghost = 0; + std::size_t n_reference = 0; + std::size_t n_contact_ghost = 0; std::size_t size = 0; // x, u, f (vector matching system dimension). @@ -408,7 +410,7 @@ class Particles auto y = Cabana::slice<0>( _aosoa_y, "current_positions" ); auto x = sliceReferencePosition(); auto u = sliceDisplacement(); - Kokkos::RangePolicy policy( 0, n_local + n_ghost ); + Kokkos::RangePolicy policy( 0, size ); auto sum_x_u = KOKKOS_LAMBDA( const std::size_t pid ) { for ( int d = 0; d < 3; d++ ) @@ -419,19 +421,22 @@ class Particles //_timer.stop(); } - void resize( int new_local, int new_ghost ) + void resize( int new_local, int new_ghost, int new_contact_ghost = 0 ) { _timer.start(); n_local = new_local; n_ghost = new_ghost; - - _plist_x.aosoa().resize( new_local + new_ghost ); - _aosoa_u.resize( new_local + new_ghost ); - _aosoa_y.resize( new_local + new_ghost ); - _aosoa_vol.resize( new_local + new_ghost ); - _plist_f.aosoa().resize( new_local ); - _aosoa_other.resize( new_local ); - _aosoa_nofail.resize( new_local + new_ghost ); + n_reference = n_local + n_ghost; + n_contact_ghost = new_contact_ghost; + int n_all = n_reference + n_contact_ghost; + + _plist_x.aosoa().resize( n_all ); + _aosoa_u.resize( n_all ); + _aosoa_y.resize( n_all ); + _aosoa_vol.resize( n_all ); + _plist_f.aosoa().resize( n_local ); + _aosoa_other.resize( n_local ); + _aosoa_nofail.resize( n_reference ); size = _plist_x.size(); _timer.stop(); }; @@ -512,9 +517,11 @@ class Particles using base_type::dim; // Per particle. + using base_type::n_contact_ghost; using base_type::n_ghost; using base_type::n_global; using base_type::n_local; + using base_type::n_reference; using base_type::size; // These are split since weighted volume only needs to be communicated once @@ -598,8 +605,8 @@ class Particles { base_type::resize( new_local, new_ghost ); _timer.start(); - _aosoa_theta.resize( new_local + new_ghost ); - _aosoa_m.resize( new_local + new_ghost ); + _aosoa_theta.resize( n_reference ); + _aosoa_m.resize( n_reference ); _timer.stop(); } @@ -673,9 +680,11 @@ class Particles using base_type::dim; // Per particle. + using base_type::n_contact_ghost; using base_type::n_ghost; using base_type::n_global; using base_type::n_local; + using base_type::n_reference; using base_type::size; // These are split since weighted volume only needs to be communicated once @@ -756,7 +765,7 @@ class Particles void resize( int new_local, int new_ghost ) { base_type::resize( new_local, new_ghost ); - _aosoa_temp.resize( new_local + new_ghost ); + _aosoa_temp.resize( n_reference ); } void output( [[maybe_unused]] const int output_step, diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 5f30dd38..18d66ba3 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -157,12 +157,6 @@ class SolverElastic // Add ghosts from other MPI ranks. comm = std::make_shared( *particles ); - std::cout << particles->n_local << " " << particles->n_ghost - << std::endl; - - if constexpr ( is_contact::value ) - contact_comm = std::make_shared( *particles ); - std::cout << particles->n_local << " " << particles->n_ghost << std::endl; @@ -200,6 +194,14 @@ class SolverElastic inputs["half_neigh"], force->_neigh_list, force_model ); } + // Purposely delay creating initial contact ghosts until after reference + // neighbor list. + if constexpr ( is_contact::value ) + contact_comm = std::make_shared( *particles ); + + std::cout << particles->n_local << " " << particles->n_ghost << " " + << particles->n_contact_ghost << std::endl; + print = print_rank(); if ( print ) { From 7018ce65f18df5aadeb32cff56f54743e826f845 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Tue, 3 Dec 2024 15:32:30 -0500 Subject: [PATCH 15/23] fixup: contact halo width --- src/CabanaPD_Comm.hpp | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/src/CabanaPD_Comm.hpp b/src/CabanaPD_Comm.hpp index 22c2686c..22d70767 100644 --- a/src/CabanaPD_Comm.hpp +++ b/src/CabanaPD_Comm.hpp @@ -400,17 +400,14 @@ class Comm std::shared_ptr gather_current; std::shared_ptr halo; HaloIds halo_ids; - int max_export; - int halo_width; // Note this initial guess is small because this is often used for very // short range interactions. - Comm( ParticleType& particles, int max_export_guess = 10 ) + Comm( ParticleType& particles, int halo_width = 1, + int max_export_guess = 10 ) : halo_ids( HaloIds( *( particles.local_grid ), particles.sliceCurrentPosition(), - particles.local_grid->haloCellWidth(), max_export_guess ) ) - , max_export( max_export_guess ) - , halo_width( particles.local_grid->haloCellWidth() ) + halo_width, max_export_guess ) ) { auto topology = Cabana::Grid::getTopology( *particles.local_grid ); halo = std::make_shared( From e311f2f65968fb302545319197a8573b298bd598 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 4 Dec 2024 10:24:17 -0500 Subject: [PATCH 16/23] Make contact optional for fragmentation problem --- examples/mechanics/fragmenting_cylinder.cpp | 32 +++++++++++++------ .../inputs/fragmenting_cylinder.json | 1 + 2 files changed, 24 insertions(+), 9 deletions(-) diff --git a/examples/mechanics/fragmenting_cylinder.cpp b/examples/mechanics/fragmenting_cylinder.cpp index 64c6dcae..94dfaa28 100644 --- a/examples/mechanics/fragmenting_cylinder.cpp +++ b/examples/mechanics/fragmenting_cylinder.cpp @@ -85,7 +85,8 @@ void fragmentingCylinderExample( const std::string filename ) { double rsq = ( x[0] - x_center ) * ( x[0] - x_center ) + ( x[1] - y_center ) * ( x[1] - y_center ); - if ( rsq < Rin * Rin || rsq > Rout * Rout || x[2] > 0.05 || x[2] < -0.05) + if ( rsq < Rin * Rin || rsq > Rout * Rout || x[2] > 0.05 || + x[2] < -0.05 ) return false; return true; }; @@ -136,17 +137,30 @@ void fragmentingCylinderExample( const std::string filename ) }; particles->updateParticles( exec_space{}, init_functor ); - double r_c = inputs["contact_horizon_factor"]; - r_c *= dx[0]; - CabanaPD::NormalRepulsionModel contact_model( delta, r_c, K ); + // ==================================================== + // Simulation run with contact physics + // ==================================================== + if ( inputs["use_contact"] ) + { + double r_c = inputs["contact_horizon_factor"]; + r_c *= dx[0]; + CabanaPD::NormalRepulsionModel contact_model( delta, r_c, K ); + auto cabana_pd = CabanaPD::createSolverFracture( + inputs, particles, force_model, contact_model ); + cabana_pd->init(); + cabana_pd->run(); + } // ==================================================== - // Simulation run + // Simulation run without contact // ==================================================== - auto cabana_pd = CabanaPD::createSolverFracture( - inputs, particles, force_model, contact_model ); - cabana_pd->init(); - cabana_pd->run(); + else + { + auto cabana_pd = CabanaPD::createSolverFracture( + inputs, particles, force_model ); + cabana_pd->init(); + cabana_pd->run(); + } } // Initialize MPI+Kokkos. diff --git a/examples/mechanics/inputs/fragmenting_cylinder.json b/examples/mechanics/inputs/fragmenting_cylinder.json index b33b28f5..cc1deba3 100644 --- a/examples/mechanics/inputs/fragmenting_cylinder.json +++ b/examples/mechanics/inputs/fragmenting_cylinder.json @@ -8,6 +8,7 @@ "shear_modulus" : {"value": 78e+9, "unit": "Pa"}, "critical_stretch" : {"value": 0.02}, "horizon" : {"value": 0.00417462, "unit": "m"}, + "use_contact" : {"value": true}, "contact_horizon_factor" : {"value": 0.9}, "cylinder_outer_radius" : {"value": 0.025, "unit": "m"}, "cylinder_inner_radius" : {"value": 0.02, "unit": "m"}, From 0b4c4887bfeed88068cf8efed903dcfe851386fa Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 4 Dec 2024 10:25:50 -0500 Subject: [PATCH 17/23] Remove impact problem --- examples/mechanics/disk_impact.cpp | 97 ---------------------- examples/mechanics/inputs/disk_impact.json | 17 ---- 2 files changed, 114 deletions(-) delete mode 100644 examples/mechanics/disk_impact.cpp delete mode 100644 examples/mechanics/inputs/disk_impact.json diff --git a/examples/mechanics/disk_impact.cpp b/examples/mechanics/disk_impact.cpp deleted file mode 100644 index bef00871..00000000 --- a/examples/mechanics/disk_impact.cpp +++ /dev/null @@ -1,97 +0,0 @@ -#include -#include -#include - -#include "mpi.h" - -#include - -#include - -// Simulate a sphere impacting a thin cylindrical plate. -void diskImpactExample( const std::string filename ) -{ - using exec_space = Kokkos::DefaultExecutionSpace; - using memory_space = typename exec_space::memory_space; - - CabanaPD::Inputs inputs( filename ); - - double K = inputs["bulk_modulus"]; - double G0 = inputs["fracture_energy"]; - double delta = inputs["horizon"]; - delta += 1e-10; - // Choose force model type. - using model_type = CabanaPD::ForceModel; - model_type force_model( delta, K, G0 ); - - // Create particles from mesh. - auto particles = CabanaPD::createParticles( - exec_space(), inputs ); - - double system_size = inputs["system_size"][0]; - auto create_func = KOKKOS_LAMBDA( const int, const double px[3] ) - { - auto radius = system_size / 2.0 - 1e-10; - auto r2 = px[0] * px[0] + px[1] * px[1]; - if ( r2 > radius * radius || px[2] < 0.0 ) - return false; - return true; - }; - particles->createParticles( exec_space{}, create_func ); - - // Define particle initialization. - auto v = particles->sliceVelocity(); - auto rho = particles->sliceDensity(); - double rho0 = inputs["density"]; - 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 ); - - double r_c = inputs["contact_horizon_factor"]; - r_c *= particles->dx[0]; - CabanaPD::NormalRepulsionModel contact_model( delta, r_c, K ); - - auto cabana_pd = CabanaPD::createSolverFracture( - inputs, particles, force_model, contact_model ); - - double impact_r = inputs["impactor_radius"]; - double impact_v = inputs["impactor_velocity"]; - double impact_x = 0.0; - double impact_y = 0.0; - auto x = particles->sliceReferencePosition(); - auto f = particles->sliceForce(); - auto body_func = KOKKOS_LAMBDA( const int p, const double t ) - { - auto z = t * impact_v + impact_r; - double r = sqrt( ( x( p, 0 ) - impact_x ) * ( x( p, 0 ) - impact_x ) + - ( x( p, 1 ) - impact_y ) * ( x( p, 1 ) - impact_y ) + - ( x( p, 2 ) - z ) * ( x( p, 2 ) - z ) ); - if ( r < impact_r ) - { - double fmag = 1.0e17 * ( r - impact_r ) * ( r - impact_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; - } - }; - auto body = CabanaPD::createBodyTerm( body_func, true ); - - cabana_pd->init( body ); - cabana_pd->run( body ); -} - -int main( int argc, char* argv[] ) -{ - MPI_Init( &argc, &argv ); - - Kokkos::initialize( argc, argv ); - - diskImpactExample( argv[1] ); - - Kokkos::finalize(); - MPI_Finalize(); -} diff --git a/examples/mechanics/inputs/disk_impact.json b/examples/mechanics/inputs/disk_impact.json deleted file mode 100644 index dc2305e7..00000000 --- a/examples/mechanics/inputs/disk_impact.json +++ /dev/null @@ -1,17 +0,0 @@ -{ - "dx" : {"value": [0.0005, 0.0005, 0.0005], "note": "This is used to enable a larger box in Z."}, - "low_corner" : {"value": [-0.03725, -0.03725, -10.0e-2], "unit": "m"}, - "high_corner" : {"value": [0.03725, 0.03725, 2.5e-3], "unit": "m"}, - "density" : {"value": 2200.0, "unit": "kg/m^3"}, - "bulk_modulus" : {"value": 14.9e+9, "unit": "Pa"}, - "fracture_energy" : {"value": 10.0575, "unit": "J/m^2"}, - "horizon" : {"value": 0.0015, "unit": "m"}, - "contact_horizon_factor" : {"value": 0.9}, - "final_time" : {"value": 2.0e-4, "unit": "s"}, - "timestep" : {"value": 1.0e-7, "unit": "s"}, - "timestep_safety_factor" : {"value": 0.85}, - "output_frequency" : {"value": 25}, - "output_reference" : {"value": false}, - "impactor_radius" : {"value": 5e-3, "unit": "m"}, - "impactor_velocity" : {"value": -100.0, "unit": "m/s"} -} From dd7a587ace3872cc78fdbc944ecfff9d9ded445a Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 4 Dec 2024 10:45:37 -0500 Subject: [PATCH 18/23] Revert contact comm --- src/CabanaPD_Comm.hpp | 89 +++++++------------------------------- src/CabanaPD_Particles.hpp | 43 +++++++----------- src/CabanaPD_Solver.hpp | 27 ------------ src/CabanaPD_Types.hpp | 3 -- 4 files changed, 32 insertions(+), 130 deletions(-) diff --git a/src/CabanaPD_Comm.hpp b/src/CabanaPD_Comm.hpp index 22d70767..6ab1508c 100644 --- a/src/CabanaPD_Comm.hpp +++ b/src/CabanaPD_Comm.hpp @@ -195,9 +195,6 @@ struct HaloIds Kokkos::parallel_for( "CabanaPD::Comm::GhostSearch", policy, ghost_search ); Kokkos::fence(); - - // Rebuild if needed. - rebuild( positions ); } template @@ -246,7 +243,6 @@ class Comm int max_export; using memory_space = typename ParticleType::memory_space; - using local_grid_type = typename ParticleType::local_grid_type; using halo_type = Cabana::Halo; using gather_u_type = Cabana::Gather; @@ -266,11 +262,12 @@ class Comm auto halo_width = local_grid->haloCellWidth(); auto topology = Cabana::Grid::getTopology( *local_grid ); - // Determine which particles should be ghosted, reallocating and - // recounting if needed. + // Determine which particles need to be ghosted to neighbors. // FIXME: set halo width based on cutoff distance. - HaloIds halo_ids( - *local_grid, positions, halo_width, max_export ); + auto halo_ids = + createHaloIds( *local_grid, positions, halo_width, max_export ); + // Rebuild if needed. + halo_ids.rebuild( positions ); // Create the Cabana Halo. halo = std::make_shared( local_grid->globalGrid().comm(), @@ -291,6 +288,17 @@ class Comm } ~Comm() {} + // Determine which particles should be ghosted, reallocating and recounting + // if needed. + template + auto createHaloIds( const LocalGridType& local_grid, + const PositionSliceType& positions, + const int min_halo_width, const int max_export ) + { + return HaloIds( + local_grid, positions, min_halo_width, max_export ); + } + // We assume here that the particle count has not changed and no resize // is necessary. void gatherDisplacement() @@ -385,71 +393,6 @@ class Comm void gatherTemperature() { gather_temp->apply(); } }; -// Does not inherit because it does not use the same reference halo -// communication pattern. -template -class Comm -{ - public: - using memory_space = typename ParticleType::memory_space; - using local_grid_type = typename ParticleType::local_grid_type; - using halo_type = Cabana::Halo; - - using gather_current_type = - Cabana::Gather; - std::shared_ptr gather_current; - std::shared_ptr halo; - HaloIds halo_ids; - - // Note this initial guess is small because this is often used for very - // short range interactions. - Comm( ParticleType& particles, int halo_width = 1, - int max_export_guess = 10 ) - : halo_ids( HaloIds( - *( particles.local_grid ), particles.sliceCurrentPosition(), - halo_width, max_export_guess ) ) - { - auto topology = Cabana::Grid::getTopology( *particles.local_grid ); - halo = std::make_shared( - particles.local_grid->globalGrid().comm(), particles.n_reference, - halo_ids._ids, halo_ids._destinations, topology ); - std::cout << "halo " << halo->numLocal() << " " << halo->numGhost() - << std::endl; - // We use n_ghost here as the "local" halo count because these current - // frame ghosts are built on top of the existing, static, reference - // frame ghosts. - particles.resize( particles.n_local, particles.n_ghost, - halo->numGhost() ); - - gather_current = - std::make_shared( *halo, particles._aosoa_y ); - } - - // This is a dynamic gather step where the steering vector needs to be - // recomputed. - void gatherCurrentPostiion( ParticleType& particles ) - { - // Get the current position. Note this is necessary to get the up to - // date current position. - auto y = particles.sliceCurrentPosition(); - // Determine which particles need to be ghosted to neighbors for the - // current positions. - halo_ids.build( y ); - - auto topology = Cabana::Grid::getTopology( *particles.local_grid ); - // FIXME: missing a build() interface - halo = std::make_shared( - particles.local_grid->globalGrid().comm(), particles.n_reference, - halo_ids._ids, halo_ids._destinations, topology ); - particles.resize( particles.n_local, particles.n_ghost, - halo->numGhost() ); - - gather_current->reserve( *halo, particles._aosoa_y ); - - gather_current->apply(); - } -}; - } // namespace CabanaPD #endif diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 24844aca..52cc36ed 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -97,8 +97,6 @@ class Particles unsigned long long int n_global = 0; std::size_t n_local = 0; std::size_t n_ghost = 0; - std::size_t n_reference = 0; - std::size_t n_contact_ghost = 0; std::size_t size = 0; // x, u, f (vector matching system dimension). @@ -139,10 +137,9 @@ class Particles // FIXME: this is for neighborlist construction. double ghost_mesh_lo[dim]; double ghost_mesh_hi[dim]; - - using local_grid_type = - Cabana::Grid::LocalGrid>; - std::shared_ptr local_grid; + std::shared_ptr< + Cabana::Grid::LocalGrid>> + local_grid; Kokkos::Array dx; int halo_width; @@ -410,7 +407,7 @@ class Particles auto y = Cabana::slice<0>( _aosoa_y, "current_positions" ); auto x = sliceReferencePosition(); auto u = sliceDisplacement(); - Kokkos::RangePolicy policy( 0, size ); + Kokkos::RangePolicy policy( 0, n_local + n_ghost ); auto sum_x_u = KOKKOS_LAMBDA( const std::size_t pid ) { for ( int d = 0; d < 3; d++ ) @@ -421,22 +418,19 @@ class Particles //_timer.stop(); } - void resize( int new_local, int new_ghost, int new_contact_ghost = 0 ) + void resize( int new_local, int new_ghost ) { _timer.start(); n_local = new_local; n_ghost = new_ghost; - n_reference = n_local + n_ghost; - n_contact_ghost = new_contact_ghost; - int n_all = n_reference + n_contact_ghost; - - _plist_x.aosoa().resize( n_all ); - _aosoa_u.resize( n_all ); - _aosoa_y.resize( n_all ); - _aosoa_vol.resize( n_all ); - _plist_f.aosoa().resize( n_local ); - _aosoa_other.resize( n_local ); - _aosoa_nofail.resize( n_reference ); + + _plist_x.aosoa().resize( new_local + new_ghost ); + _aosoa_u.resize( new_local + new_ghost ); + _aosoa_y.resize( new_local + new_ghost ); + _aosoa_vol.resize( new_local + new_ghost ); + _plist_f.aosoa().resize( new_local ); + _aosoa_other.resize( new_local ); + _aosoa_nofail.resize( new_local + new_ghost ); size = _plist_x.size(); _timer.stop(); }; @@ -482,7 +476,6 @@ class Particles friend class Comm; friend class Comm; - friend class Comm; protected: aosoa_u_type _aosoa_u; @@ -517,11 +510,9 @@ class Particles using base_type::dim; // Per particle. - using base_type::n_contact_ghost; using base_type::n_ghost; using base_type::n_global; using base_type::n_local; - using base_type::n_reference; using base_type::size; // These are split since weighted volume only needs to be communicated once @@ -605,8 +596,8 @@ class Particles { base_type::resize( new_local, new_ghost ); _timer.start(); - _aosoa_theta.resize( n_reference ); - _aosoa_m.resize( n_reference ); + _aosoa_theta.resize( new_local + new_ghost ); + _aosoa_m.resize( new_local + new_ghost ); _timer.stop(); } @@ -680,11 +671,9 @@ class Particles using base_type::dim; // Per particle. - using base_type::n_contact_ghost; using base_type::n_ghost; using base_type::n_global; using base_type::n_local; - using base_type::n_reference; using base_type::size; // These are split since weighted volume only needs to be communicated once @@ -765,7 +754,7 @@ class Particles void resize( int new_local, int new_ghost ) { base_type::resize( new_local, new_ghost ); - _aosoa_temp.resize( n_reference ); + _aosoa_temp.resize( new_local + new_ghost ); } void output( [[maybe_unused]] const int output_step, diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 18d66ba3..a1db89a2 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -113,8 +113,6 @@ class SolverElastic using heat_transfer_type = HeatTransfer; using contact_type = Force; using contact_model_type = ContactModel; - using contact_comm_type = - Comm; SolverElastic( input_type _inputs, std::shared_ptr _particles, @@ -152,14 +150,9 @@ class SolverElastic dt = inputs["timestep"]; integrator = std::make_shared( dt ); - std::cout << particles->n_local << " " << particles->n_ghost - << std::endl; // Add ghosts from other MPI ranks. comm = std::make_shared( *particles ); - std::cout << particles->n_local << " " << particles->n_ghost - << std::endl; - // Update temperature ghost size if needed. if constexpr ( is_temperature_dependent< typename force_model_type::thermal_type>::value ) @@ -194,14 +187,6 @@ class SolverElastic inputs["half_neigh"], force->_neigh_list, force_model ); } - // Purposely delay creating initial contact ghosts until after reference - // neighbor list. - if constexpr ( is_contact::value ) - contact_comm = std::make_shared( *particles ); - - std::cout << particles->n_local << " " << particles->n_ghost << " " - << particles->n_contact_ghost << std::endl; - print = print_rank(); if ( print ) { @@ -477,7 +462,6 @@ class SolverElastic // Optional modules. std::shared_ptr heat_transfer; std::shared_ptr contact; - std::shared_ptr contact_comm; contact_model_type contact_model; // Output files. @@ -626,11 +610,6 @@ class SolverFracture // Update ghost particles. comm->gatherDisplacement(); - // Need to ghost current positions for DEM/contact. - if constexpr ( is_contact::value || - is_contact::value ) - contact_comm->gatherCurrentPostiion( *particles ); - // Compute internal forces. updateForce(); @@ -670,11 +649,6 @@ class SolverFracture // Update ghost particles. comm->gatherDisplacement(); - // Need to ghost current positions for DEM/contact. - if constexpr ( is_contact::value || - is_contact::value ) - contact_comm->gatherCurrentPostiion( *particles ); - // Compute internal forces. updateForce(); @@ -734,7 +708,6 @@ class SolverFracture protected: using base_type::comm; using base_type::contact; - using base_type::contact_comm; using base_type::force; using base_type::inputs; using base_type::integrator; diff --git a/src/CabanaPD_Types.hpp b/src/CabanaPD_Types.hpp index 2e2c5c98..4ca10d35 100644 --- a/src/CabanaPD_Types.hpp +++ b/src/CabanaPD_Types.hpp @@ -24,9 +24,6 @@ struct Fracture // Contact and DEM (contact without PD) tags. -struct Contact -{ -}; struct NoContact { }; From 8182615defad527e3f1ad9888f65e40572eca7d4 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 4 Dec 2024 11:01:42 -0500 Subject: [PATCH 19/23] fixup: remove --- examples/mechanics/CMakeLists.txt | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/examples/mechanics/CMakeLists.txt b/examples/mechanics/CMakeLists.txt index 19148fa3..a10b05f1 100644 --- a/examples/mechanics/CMakeLists.txt +++ b/examples/mechanics/CMakeLists.txt @@ -10,7 +10,4 @@ target_link_libraries(CrackBranching LINK_PUBLIC CabanaPD) add_executable(FragmentingCylinder fragmenting_cylinder.cpp) target_link_libraries(FragmentingCylinder LINK_PUBLIC CabanaPD) -add_executable(DiskImpact disk_impact.cpp) -target_link_libraries(DiskImpact LINK_PUBLIC CabanaPD) - -install(TARGETS ElasticWave KalthoffWinkler CrackBranching FragmentingCylinder DiskImpact DESTINATION ${CMAKE_INSTALL_BINDIR}) +install(TARGETS ElasticWave KalthoffWinkler CrackBranching FragmentingCylinder DESTINATION ${CMAKE_INSTALL_BINDIR}) From b0fd1fa646dfdd3134551bb70859702f09843b32 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 4 Dec 2024 11:02:39 -0500 Subject: [PATCH 20/23] Disable contact with MPI --- src/CabanaPD_Comm.hpp | 4 +++- src/CabanaPD_Solver.hpp | 7 +++++++ 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/src/CabanaPD_Comm.hpp b/src/CabanaPD_Comm.hpp index 6ab1508c..5ee56eb3 100644 --- a/src/CabanaPD_Comm.hpp +++ b/src/CabanaPD_Comm.hpp @@ -286,7 +286,9 @@ class Comm _init_timer.stop(); } - ~Comm() {} + + auto size() { return mpi_size; } + auto rank() { return mpi_rank; } // Determine which particles should be ghosted, reallocating and recounting // if needed. diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index a1db89a2..37469a1d 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -153,6 +153,13 @@ class SolverElastic // Add ghosts from other MPI ranks. comm = std::make_shared( *particles ); + if constexpr ( is_contact::value ) + { + if ( comm->size() > 1 ) + throw std::runtime_error( + "Contact with MPI is currently disabled." ); + } + // Update temperature ghost size if needed. if constexpr ( is_temperature_dependent< typename force_model_type::thermal_type>::value ) From 5037bdb93b51cc5618861bc460e2ed1ed1a01b2b Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Thu, 5 Dec 2024 09:20:54 -0500 Subject: [PATCH 21/23] fixup: simplify solver creation functions --- src/CabanaPD_Solver.hpp | 39 ++++++++++++++++++++++++++++++--------- 1 file changed, 30 insertions(+), 9 deletions(-) diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 37469a1d..6379734f 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -732,26 +732,47 @@ class SolverFracture // =============================================================== template + class ForceModel> auto createSolverElastic( InputsType inputs, std::shared_ptr particles, - ForceModel model, ContactModelType... contact_model ) + ForceModel model ) +{ + return std::make_shared< + SolverElastic>( + inputs, particles, model ); +} + +template +auto createSolverElastic( InputsType inputs, + std::shared_ptr particles, + ForceModel model, ContactModelType contact_model ) { return std::make_shared>( - inputs, particles, model, contact_model... ); + ForceModel, ContactModelType>>( + inputs, particles, model, contact_model ); } template + class ForceModel> auto createSolverFracture( InputsType inputs, std::shared_ptr particles, - ForceModel model, ContactModelType... contact_model ) + ForceModel model ) { return std::make_shared< - SolverFracture>( inputs, particles, model, - contact_model... ); + SolverFracture>( + inputs, particles, model ); +} + +template +auto createSolverFracture( InputsType inputs, + std::shared_ptr particles, + ForceModel model, ContactModelType contact_model ) +{ + return std::make_shared>( + inputs, particles, model, contact_model ); } } // namespace CabanaPD From af9b13613d1f828421b51bd3ff715c9cb6a1279f Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Thu, 5 Dec 2024 09:21:02 -0500 Subject: [PATCH 22/23] fixup: rename contact models file --- src/CabanaPD.hpp | 2 +- ...anaPD_ForceModels_Contact.hpp => CabanaPD_ContactModels.hpp} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename src/force/{CabanaPD_ForceModels_Contact.hpp => CabanaPD_ContactModels.hpp} (100%) diff --git a/src/CabanaPD.hpp b/src/CabanaPD.hpp index 4acf096d..f1014885 100644 --- a/src/CabanaPD.hpp +++ b/src/CabanaPD.hpp @@ -29,7 +29,7 @@ #include #include -#include +#include #include #include #include diff --git a/src/force/CabanaPD_ForceModels_Contact.hpp b/src/force/CabanaPD_ContactModels.hpp similarity index 100% rename from src/force/CabanaPD_ForceModels_Contact.hpp rename to src/force/CabanaPD_ContactModels.hpp From 71cf21ea34fe15452b5c08c88b853c06215bc473 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Thu, 5 Dec 2024 11:33:41 -0500 Subject: [PATCH 23/23] fixup: rename template model types --- src/CabanaPD_Solver.hpp | 50 +++++++++++++++++++++-------------------- 1 file changed, 26 insertions(+), 24 deletions(-) diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 6379734f..d8be3db4 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -92,7 +92,7 @@ class SolverBase }; template + class ForceModelType, class ContactModelType = NoContact> class SolverElastic { public: @@ -102,7 +102,7 @@ class SolverElastic // Core module types - required for all problems. using particle_type = ParticleType; using integrator_type = Integrator; - using force_model_type = ForceModel; + using force_model_type = ForceModelType; using force_type = Force; using comm_type = Comm; @@ -111,8 +111,8 @@ class SolverElastic // Optional module types. using heat_transfer_type = HeatTransfer; - using contact_type = Force; - using contact_model_type = ContactModel; + using contact_type = Force; + using contact_model_type = ContactModelType; SolverElastic( input_type _inputs, std::shared_ptr _particles, @@ -485,26 +485,26 @@ class SolverElastic }; template + class ForceModelType, class ContactModelType = NoContact> class SolverFracture - : public SolverElastic + : public SolverElastic { public: using base_type = SolverElastic; + ForceModelType, ContactModelType>; using exec_space = typename base_type::exec_space; using memory_space = typename base_type::memory_space; using particle_type = typename base_type::particle_type; using integrator_type = typename base_type::integrator_type; using comm_type = typename base_type::comm_type; - using force_model_type = ForceModel; + using force_model_type = ForceModelType; using force_type = typename base_type::force_type; using neigh_iter_tag = Cabana::SerialOpTag; using input_type = typename base_type::input_type; - using contact_model_type = ContactModel; + using contact_model_type = ContactModelType; SolverFracture( input_type _inputs, std::shared_ptr _particles, @@ -732,47 +732,49 @@ class SolverFracture // =============================================================== template + class ForceModelType> auto createSolverElastic( InputsType inputs, std::shared_ptr particles, - ForceModel model ) + ForceModelType model ) { return std::make_shared< - SolverElastic>( + SolverElastic>( inputs, particles, model ); } template + class ForceModelType, class ContactModelType> auto createSolverElastic( InputsType inputs, std::shared_ptr particles, - ForceModel model, ContactModelType contact_model ) + ForceModelType model, ContactModelType contact_model ) { return std::make_shared>( + ForceModelType, ContactModelType>>( inputs, particles, model, contact_model ); } template + class ForceModelType> auto createSolverFracture( InputsType inputs, std::shared_ptr particles, - ForceModel model ) + ForceModelType model ) { return std::make_shared< - SolverFracture>( + SolverFracture>( inputs, particles, model ); } template + class ForceModelType, class ContactModelType> auto createSolverFracture( InputsType inputs, std::shared_ptr particles, - ForceModel model, ContactModelType contact_model ) + ForceModelType model, + ContactModelType contact_model ) { - return std::make_shared>( - inputs, particles, model, contact_model ); + return std::make_shared< + SolverFracture>( inputs, particles, model, + contact_model ); } } // namespace CabanaPD