From 7f41184ba4b83a2c5f9b88905082de163e93cc28 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Thu, 11 Jul 2024 10:46:11 -0400 Subject: [PATCH 1/7] Remove unused default constructors and set_param --- src/CabanaPD_ForceModels.hpp | 9 ----- src/force/CabanaPD_ForceModels_LPS.hpp | 22 ++--------- src/force/CabanaPD_ForceModels_PMB.hpp | 53 ++------------------------ 3 files changed, 7 insertions(+), 77 deletions(-) diff --git a/src/CabanaPD_ForceModels.hpp b/src/CabanaPD_ForceModels.hpp index 6685de1c..a5c43790 100644 --- a/src/CabanaPD_ForceModels.hpp +++ b/src/CabanaPD_ForceModels.hpp @@ -20,12 +20,9 @@ struct BaseForceModel { double delta; - BaseForceModel(){}; BaseForceModel( const double _delta ) : delta( _delta ){}; - BaseForceModel( const BaseForceModel& model ) { delta = model.delta; } - // No-op for temperature. KOKKOS_INLINE_FUNCTION void thermalStretch( double&, const int, const int ) const {} @@ -54,12 +51,6 @@ struct BaseTemperatureModel temperature = model.temperature; } - void set_param( const double _alpha, const double _temp0 ) - { - alpha = _alpha; - temp0 = _temp0; - } - void update( const TemperatureType _temp ) { temperature = _temp; } // Update stretch with temperature effects. diff --git a/src/force/CabanaPD_ForceModels_LPS.hpp b/src/force/CabanaPD_ForceModels_LPS.hpp index 8f2999fc..f009e01b 100644 --- a/src/force/CabanaPD_ForceModels_LPS.hpp +++ b/src/force/CabanaPD_ForceModels_LPS.hpp @@ -34,21 +34,13 @@ struct ForceModel : public BaseForceModel double theta_coeff; double s_coeff; - ForceModel(){}; ForceModel( const double _delta, const double _K, const double _G, const int _influence = 0 ) : base_type( _delta ) , influence_type( _influence ) + , K( _K ) + , G( _G ) { - set_param( _delta, _K, _G ); - } - - void set_param( const double _delta, const double _K, const double _G ) - { - delta = _delta; - K = _K; - G = _G; - theta_coeff = 3.0 * K - 5.0 * G; s_coeff = 15.0 * G; } @@ -80,19 +72,11 @@ struct ForceModel : public ForceModel double s0; double bond_break_coeff; - ForceModel() {} ForceModel( const double _delta, const double _K, const double _G, const double _G0, const int _influence = 0 ) : base_type( _delta, _K, _G, _influence ) + , G0( _G0 ) { - set_param( _delta, _K, _G, _G0 ); - } - - void set_param( const double _delta, const double _K, const double _G, - const double _G0 ) - { - base_type::set_param( _delta, _K, _G ); - G0 = _G0; if ( influence_type == 1 ) { s0 = sqrt( 5.0 * G0 / 9.0 / K / delta ); // 1/xi diff --git a/src/force/CabanaPD_ForceModels_PMB.hpp b/src/force/CabanaPD_ForceModels_PMB.hpp index c80acd63..15deca9e 100644 --- a/src/force/CabanaPD_ForceModels_PMB.hpp +++ b/src/force/CabanaPD_ForceModels_PMB.hpp @@ -31,24 +31,10 @@ struct ForceModel : public BaseForceModel double c; double K; - ForceModel(){}; - ForceModel( const double delta, const double K ) + ForceModel( const double delta, const double _K ) : base_type( delta ) + , K( _K ) { - set_param( delta, K ); - } - - ForceModel( const ForceModel& model ) - : base_type( model ) - { - c = model.c; - K = model.K; - } - - void set_param( const double _delta, const double _K ) - { - delta = _delta; - K = _K; c = 18.0 * K / ( pi * delta * delta * delta * delta ); } }; @@ -69,25 +55,10 @@ struct ForceModel double s0; double bond_break_coeff; - ForceModel() {} - ForceModel( const double delta, const double K, const double G0 ) + ForceModel( const double delta, const double K, const double _G0 ) : base_type( delta, K ) + , G0( _G0 ) { - set_param( delta, K, G0 ); - } - - ForceModel( const ForceModel& model ) - : base_type( model ) - { - G0 = model.G0; - s0 = model.s0; - bond_break_coeff = model.bond_break_coeff; - } - - void set_param( const double _delta, const double _K, const double _G0 ) - { - base_type::set_param( _delta, _K ); - G0 = _G0; s0 = sqrt( 5.0 * G0 / 9.0 / K / delta ); bond_break_coeff = ( 1.0 + s0 ) * ( 1.0 + s0 ); } @@ -165,14 +136,6 @@ struct ForceModel : base_type( _delta, _K ) , base_temperature_type( _temp, _alpha, _temp0 ) { - set_param( _delta, _K, _alpha, _temp0 ); - } - - void set_param( const double _delta, const double _K, const double _alpha, - const double _temp0 ) - { - base_type::set_param( _delta, _K ); - base_temperature_type::set_param( _alpha, _temp0 ); } }; @@ -220,14 +183,6 @@ struct ForceModel : base_type( _delta, _K, _G0 ) , base_temperature_type( _temp, _alpha, _temp0 ) { - set_param( _delta, _K, _G0, _alpha, _temp0 ); - } - - void set_param( const double _delta, const double _K, const double _G0, - const double _alpha, const double _temp0 ) - { - base_type::set_param( _delta, _K, _G0 ); - base_temperature_type::set_param( _alpha, _temp0 ); } KOKKOS_INLINE_FUNCTION From faa55d8af1d005519438a31c036a0b042b7517aa Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Tue, 16 Jul 2024 22:17:03 -0400 Subject: [PATCH 2/7] Add support for multiple types/species/materials --- .../thermomechanics/thermal_deformation.cpp | 3 +- src/CabanaPD_Comm.hpp | 58 +++- src/CabanaPD_Force.hpp | 5 +- src/CabanaPD_ForceModels.hpp | 69 ++++- src/CabanaPD_Particles.hpp | 41 +-- src/CabanaPD_Solver.hpp | 39 ++- src/CabanaPD_Types.hpp | 8 + src/force/CabanaPD_ForceModels_LPS.hpp | 9 +- src/force/CabanaPD_ForceModels_PMB.hpp | 258 ++++++++++++++++-- src/force/CabanaPD_Force_LPS.hpp | 44 +-- src/force/CabanaPD_Force_PMB.hpp | 25 +- unit_test/tstComm.hpp | 2 +- 12 files changed, 468 insertions(+), 93 deletions(-) diff --git a/examples/thermomechanics/thermal_deformation.cpp b/examples/thermomechanics/thermal_deformation.cpp index 4e027a71..488cc2fd 100644 --- a/examples/thermomechanics/thermal_deformation.cpp +++ b/examples/thermomechanics/thermal_deformation.cpp @@ -82,7 +82,8 @@ void thermalDeformationExample( const std::string filename ) // Force model // ==================================================== auto force_model = CabanaPD::createForceModel( - model_type{}, CabanaPD::Elastic{}, *particles, delta, K, alpha, temp0 ); + model_type{}, CabanaPD::Elastic{}, CabanaPD::TemperatureDependent{}, + *particles, delta, K, alpha, temp0 ); // ==================================================== // Create solver diff --git a/src/CabanaPD_Comm.hpp b/src/CabanaPD_Comm.hpp index 6ab1508c..3b73cd34 100644 --- a/src/CabanaPD_Comm.hpp +++ b/src/CabanaPD_Comm.hpp @@ -230,12 +230,13 @@ struct HaloIds } }; -template +template class Comm; // FIXME: extract model from ParticleType instead. template -class Comm +class Comm { public: int mpi_size = -1; @@ -320,11 +321,12 @@ class Comm }; template -class Comm - : public Comm +class Comm + : public Comm { public: - using base_type = Comm; + using base_type = + Comm; using memory_space = typename base_type::memory_space; using halo_type = typename base_type::halo_type; using base_type::gather_u; @@ -368,12 +370,52 @@ class Comm } }; +template +class Comm + : public Comm +{ + public: + using base_type = + Comm; + using memory_space = typename base_type::memory_space; + using halo_type = typename base_type::halo_type; + using base_type::halo; + + using base_type::_init_timer; + using base_type::_timer; + + using gather_species_type = + Cabana::Gather; + std::shared_ptr gather_species; + + Comm( ParticleType& particles, int max_export_guess = 100 ) + : base_type( particles, max_export_guess ) + { + _init_timer.start(); + + gather_species = std::make_shared( + *halo, particles._aosoa_species ); + + particles.resize( halo->numLocal(), halo->numGhost() ); + _init_timer.stop(); + } + + void gatherSpecies() + { + _timer.start(); + gather_species->apply(); + _timer.stop(); + } +}; + template -class Comm - : public Comm +class Comm + : public Comm { public: - using base_type = Comm; + using base_type = + Comm; using memory_space = typename base_type::memory_space; using halo_type = typename base_type::halo_type; using base_type::halo; diff --git a/src/CabanaPD_Force.hpp b/src/CabanaPD_Force.hpp index 9ea0e040..77707ef1 100644 --- a/src/CabanaPD_Force.hpp +++ b/src/CabanaPD_Force.hpp @@ -129,8 +129,7 @@ getLinearizedDistance( const PosType& x, const PosType& u, const int i, template class Force; -template -class Force +class ForceBase { protected: bool _half_neigh; @@ -139,7 +138,7 @@ class Force Timer _energy_timer; public: - Force( const bool half_neigh ) + ForceBase( const bool half_neigh ) : _half_neigh( half_neigh ) { } diff --git a/src/CabanaPD_ForceModels.hpp b/src/CabanaPD_ForceModels.hpp index a5c43790..93afb9ca 100644 --- a/src/CabanaPD_ForceModels.hpp +++ b/src/CabanaPD_ForceModels.hpp @@ -16,13 +16,71 @@ namespace CabanaPD { -struct BaseForceModel + +template +struct BaseForceModel; + +template <> +struct BaseForceModel<> { + using species_type = SingleSpecies; double delta; + BaseForceModel(){}; BaseForceModel( const double _delta ) : delta( _delta ){}; + auto horizon( const int ) { return delta; } + auto maxHorizon() { return delta; } + + // No-op for temperature. + KOKKOS_INLINE_FUNCTION + void thermalStretch( double&, const int, const int ) const {} +}; + +template +struct BaseForceModel +{ + using species_type = MultiSpecies; + + // Only allow one memory space. + using memory_space = + typename std::tuple_element<0, std::tuple>::type; + using view_type_1d = Kokkos::View; + view_type_1d delta; + double max_delta; + std::size_t num_types; + + using view_type_2d = Kokkos::View; + + BaseForceModel( const double _delta ) + : delta( view_type_1d( "delta", 1 ) ) + , num_types( 1 ) + { + Kokkos::deep_copy( delta, _delta ); + } + + template + BaseForceModel( const ArrayType& _delta ) + : delta( view_type_1d( "delta", _delta.size() ) ) + , num_types( _delta.size() ) + { + max_delta = 0; + auto init_func = KOKKOS_LAMBDA( const int i, double& max ) + { + delta( i ) = _delta[i]; + if ( delta( i ) > max ) + max = delta( i ); + }; + using exec_space = typename memory_space::execution_space; + Kokkos::RangePolicy policy( 0, num_types ); + Kokkos::parallel_reduce( "CabanaPD::Model::Init", policy, init_func, + max_delta ); + } + + auto horizon( const int i ) { return delta( i ); } + auto maxHorizon() { return max_delta; } + // No-op for temperature. KOKKOS_INLINE_FUNCTION void thermalStretch( double&, const int, const int ) const {} @@ -31,6 +89,8 @@ struct BaseForceModel template struct BaseTemperatureModel { + using species_type = SingleSpecies; + using memory_space = typename TemperatureType::memory_space; double alpha; double temp0; @@ -44,13 +104,6 @@ struct BaseTemperatureModel , temp0( _temp0 ) , temperature( _temp ){}; - BaseTemperatureModel( const BaseTemperatureModel& model ) - { - alpha = model.alpha; - temp0 = model.temp0; - temperature = model.temperature; - } - void update( const TemperatureType _temp ) { temperature = _temp; } // Update stretch with temperature effects. diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index c9e9816d..d545127c 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -101,12 +101,11 @@ class Particles using vector_type = Cabana::MemberTypes; // volume, dilatation, weighted_volume. using scalar_type = Cabana::MemberTypes; - // no-fail. + // no-fail, type. using int_type = Cabana::MemberTypes; - // v, W, rho, damage, type. + // v, W, rho, damage. using other_types = - Cabana::MemberTypes; - // Potentially needed later: body force (b), ID. + Cabana::MemberTypes; // FIXME: add vector length. // FIXME: enable variable aosoa. @@ -114,6 +113,7 @@ class Particles using aosoa_y_type = Cabana::AoSoA; using aosoa_vol_type = Cabana::AoSoA; using aosoa_nofail_type = Cabana::AoSoA; + using aosoa_species_type = Cabana::AoSoA; using aosoa_other_type = Cabana::AoSoA; // Using grid here for the particle init. using plist_x_type = @@ -356,8 +356,11 @@ class Particles { return Cabana::slice<0>( _aosoa_vol, "volume" ); } - auto sliceType() { return Cabana::slice<4>( _aosoa_other, "type" ); } - auto sliceType() const { return Cabana::slice<4>( _aosoa_other, "type" ); } + auto sliceType() { return Cabana::slice<0>( _aosoa_species, "type" ); } + auto sliceType() const + { + return Cabana::slice<0>( _aosoa_species, "type" ); + } auto sliceStrainEnergy() { return Cabana::slice<1>( _aosoa_other, "strain_energy" ); @@ -427,6 +430,7 @@ class Particles _plist_f.aosoa().resize( new_local ); _aosoa_other.resize( new_local ); _aosoa_nofail.resize( new_local + new_ghost ); + _aosoa_species.resize( new_local + new_ghost ); size = _plist_x.size(); _timer.stop(); }; @@ -457,7 +461,7 @@ class Particles "particles", local_grid->globalGrid(), output_step, output_time, 0, n_local, getPosition( use_reference ), sliceStrainEnergy(), sliceForce(), sliceDisplacement(), sliceVelocity(), - sliceDamage() ); + sliceDamage(), sliceType() ); #else log( std::cout, "No particle output enabled." ); #endif @@ -470,14 +474,17 @@ class Particles auto timeOutput() { return _output_timer.time(); }; auto time() { return _timer.time(); }; - friend class Comm; - friend class Comm; + friend class Comm; + friend class Comm; + friend class Comm; + friend class Comm; protected: aosoa_u_type _aosoa_u; aosoa_y_type _aosoa_y; aosoa_vol_type _aosoa_vol; aosoa_nofail_type _aosoa_nofail; + aosoa_species_type _aosoa_species; aosoa_other_type _aosoa_other; plist_x_type _plist_x; @@ -619,7 +626,7 @@ class Particles base_type::sliceStrainEnergy(), base_type::sliceForce(), base_type::sliceDisplacement(), base_type::sliceVelocity(), base_type::sliceDamage(), sliceWeightedVolume(), - sliceDilatation() ); + sliceDilatation(), base_type::sliceType() ); #else log( std::cout, "No particle output enabled." ); #endif @@ -628,8 +635,10 @@ class Particles _output_timer.stop(); } - friend class Comm; - friend class Comm; + friend class Comm; + friend class Comm; + friend class Comm; + friend class Comm; protected: void init_lps() @@ -761,10 +770,10 @@ class Particles #endif } - friend class Comm; - friend class Comm; - friend class Comm; - friend class Comm; + friend class Comm; + friend class Comm; + friend class Comm; + friend class Comm; protected: void init_temp() diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 04329965..815a6152 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -103,6 +103,7 @@ class SolverElastic using force_model_type = ForceModel; using force_type = Force; using comm_type = Comm; using neighbor_type = Cabana::VerletList( *particles ); - // Update temperature ghost size if needed. + // Update optional property ghost sizes if needed. + if constexpr ( std::is_same::value ) + force_model.update( particles->sliceType() ); if constexpr ( std::is_same::value ) force_model.update( particles->sliceTemperature() ); @@ -146,8 +150,8 @@ class SolverElastic particles->ghost_mesh_hi[2] }; auto x = particles->sliceReferencePosition(); neighbors = std::make_shared( x, 0, particles->n_local, - force_model.delta, 1.0, - mesh_min, mesh_max ); + force_model.maxHorizon(), + 1.0, mesh_min, mesh_max ); _neighbor_timer.stop(); _init_timer.start(); @@ -229,7 +233,10 @@ class SolverElastic if ( !boundary_condition.forceUpdate() ) boundary_condition.apply( exec_space(), *particles, 0.0 ); - // Communicate temperature. + // Communicate optional properties. + if constexpr ( std::is_same::value ) + comm->gatherSpecies(); if constexpr ( std::is_same::value ) comm->gatherTemperature(); @@ -269,6 +276,11 @@ class SolverElastic // Update ghost particles. comm->gatherDisplacement(); + // Communicate optional type. + if constexpr ( std::is_same::value ) + comm->gatherSpecies(); + // Compute internal forces. updateForce(); @@ -301,6 +313,11 @@ class SolverElastic // Update ghost particles. comm->gatherDisplacement(); + // Communicate optional properties. + if constexpr ( std::is_same::value ) + comm->gatherSpecies(); + // Compute internal forces. updateForce(); @@ -522,7 +539,10 @@ class SolverFracture if ( !boundary_condition.forceUpdate() ) boundary_condition.apply( exec_space(), *particles, 0.0 ); - // Communicate temperature. + // Communicate optional properties. + if constexpr ( std::is_same::value ) + comm->gatherSpecies(); if constexpr ( std::is_same::value ) comm->gatherTemperature(); @@ -562,6 +582,11 @@ class SolverFracture // Update ghost particles. comm->gatherDisplacement(); + // Communicate optional type. + if constexpr ( std::is_same::value ) + comm->gatherSpecies(); + // Compute internal forces. updateForce(); @@ -597,6 +622,10 @@ class SolverFracture // Update ghost particles. comm->gatherDisplacement(); + // Communicate optional properties. + if constexpr ( std::is_same::value ) + comm->gatherSpecies(); // Compute internal forces. updateForce(); diff --git a/src/CabanaPD_Types.hpp b/src/CabanaPD_Types.hpp index 62076735..ac5d097a 100644 --- a/src/CabanaPD_Types.hpp +++ b/src/CabanaPD_Types.hpp @@ -20,6 +20,14 @@ struct Elastic struct Fracture { }; + +struct SingleSpecies +{ +}; +struct MultiSpecies +{ +}; + struct TemperatureDependent { }; diff --git a/src/force/CabanaPD_ForceModels_LPS.hpp b/src/force/CabanaPD_ForceModels_LPS.hpp index f009e01b..14bd9be8 100644 --- a/src/force/CabanaPD_ForceModels_LPS.hpp +++ b/src/force/CabanaPD_ForceModels_LPS.hpp @@ -17,10 +17,12 @@ namespace CabanaPD { + template <> -struct ForceModel : public BaseForceModel +struct ForceModel : public BaseForceModel<> { - using base_type = BaseForceModel; + using base_type = BaseForceModel<>; + using species_type = SingleSpecies; using base_model = LPS; using fracture_type = Elastic; using thermal_type = TemperatureIndependent; @@ -58,6 +60,7 @@ template <> struct ForceModel : public ForceModel { using base_type = ForceModel; + using species_type = SingleSpecies; using base_model = typename base_type::base_model; using fracture_type = Fracture; using thermal_type = base_type::thermal_type; @@ -93,6 +96,7 @@ template <> struct ForceModel : public ForceModel { using base_type = ForceModel; + using species_type = SingleSpecies; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; using thermal_type = base_type::thermal_type; @@ -111,6 +115,7 @@ template <> struct ForceModel : public ForceModel { using base_type = ForceModel; + using species_type = SingleSpecies; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; using thermal_type = base_type::thermal_type; diff --git a/src/force/CabanaPD_ForceModels_PMB.hpp b/src/force/CabanaPD_ForceModels_PMB.hpp index 15deca9e..915fa9d0 100644 --- a/src/force/CabanaPD_ForceModels_PMB.hpp +++ b/src/force/CabanaPD_ForceModels_PMB.hpp @@ -18,10 +18,13 @@ namespace CabanaPD { + template <> -struct ForceModel : public BaseForceModel +struct ForceModel + : public BaseForceModel<> { - using base_type = BaseForceModel; + using base_type = BaseForceModel<>; + using species_type = typename base_type::species_type; using base_model = PMB; using fracture_type = Elastic; using thermal_type = TemperatureIndependent; @@ -31,14 +34,111 @@ struct ForceModel : public BaseForceModel double c; double K; - ForceModel( const double delta, const double _K ) + ForceModel(){}; + ForceModel( const double delta, const double K ) + : base_type( delta ) + { + setParameters( delta, K ); + } + + void setParameters( const double _delta, const double _K ) + { + delta = _delta; + K = _K; + c = micromodulus(); + } + + KOKKOS_INLINE_FUNCTION + double micromodulus() + { + return 18.0 * K / ( pi * delta * delta * delta * delta ); + } + + KOKKOS_INLINE_FUNCTION + auto micromodulus( const int, const int ) const { return c; } +}; + +template +struct ForceModel + : public BaseForceModel +{ + using memory_space = typename ParticleType::memory_space; + using base_type = BaseForceModel; + using species_type = typename base_type::species_type; + using base_model = PMB; + using fracture_type = Elastic; + using thermal_type = TemperatureIndependent; + + using base_type::delta; + using base_type::num_types; + using view_type_1d = typename base_type::view_type_1d; + using view_type_2d = typename base_type::view_type_2d; + view_type_2d c; + view_type_1d K; + ParticleType type; + + template + ForceModel( const ArrayType& delta, const ArrayType& _K, + ParticleType _type ) : base_type( delta ) - , K( _K ) + , c( view_type_2d( "micromodulus", num_types, num_types ) ) + , K( view_type_1d( "bulk_modulus", num_types ) ) + , type( _type ) + { + setParameters( _K ); + } + + template + void setParameters( const ArrayType& _K ) + { + // Initialize self interaction parameters. + auto init_self_func = KOKKOS_LAMBDA( const int i ) + { + K( i ) = _K[i]; + c( i, i ) = micromodulus( i ); + }; + using exec_space = typename memory_space::execution_space; + Kokkos::RangePolicy policy( 0, num_types ); + Kokkos::parallel_for( "CabanaPD::Model::Init", policy, init_self_func ); + Kokkos::fence(); + + // Initialize cross-terms. + auto init_cross_func = KOKKOS_LAMBDA( const int i ) + { + for ( std::size_t j = i; j < num_types; j++ ) + c( i, j ) = ( micromodulus( i ) + micromodulus( j ) ) / 2.0; + }; + Kokkos::parallel_for( "CabanaPD::Model::Init", policy, + init_cross_func ); + } + + KOKKOS_INLINE_FUNCTION + auto micromodulus( const int i ) const + { + auto d = delta( i ); + return 18.0 * K( i ) / ( pi * d * d * d * d ); + } + + KOKKOS_INLINE_FUNCTION + auto micromodulus( const int i, const int j ) const { - c = 18.0 * K / ( pi * delta * delta * delta * delta ); + return c( type( i ), type( j ) ); } + + void update( const ParticleType _type ) { type = _type; } }; +template +auto createForceModel( PMB, Elastic, TemperatureIndependent, + ParticleType particles, const double delta, + const double K ) +{ + auto type = particles.sliceType(); + using type_type = decltype( type ); + return ForceModel( + delta, K, type ); +} + template <> struct ForceModel : public ForceModel @@ -51,14 +151,21 @@ struct ForceModel using base_type::c; using base_type::delta; using base_type::K; + double G0; double s0; double bond_break_coeff; - ForceModel( const double delta, const double K, const double _G0 ) + ForceModel() {} + ForceModel( const double delta, const double K, const double G0 ) : base_type( delta, K ) - , G0( _G0 ) { + setParameters( G0 ); + } + + void setParameters( const double _G0 ) + { + G0 = _G0; s0 = sqrt( 5.0 * G0 / 9.0 / K / delta ); bond_break_coeff = ( 1.0 + s0 ) * ( 1.0 + s0 ); } @@ -71,14 +178,122 @@ struct ForceModel } }; -template <> -struct ForceModel - : public ForceModel +template +struct ForceModel + : public ForceModel { - using base_type = ForceModel; + using base_type = + ForceModel; + using memory_space = typename base_type::memory_space; + using base_model = typename base_type::base_model; + using species_type = typename base_type::species_type; + using fracture_type = Fracture; + using thermal_type = typename base_type::thermal_type; + + using base_type::c; + using base_type::delta; + using base_type::K; + using base_type::num_types; + using base_type::type; + + using view_type_1d = typename base_type::view_type_1d; + using view_type_2d = typename base_type::view_type_2d; + view_type_1d G0; + view_type_2d s0; + view_type_2d bond_break_coeff; + + template + ForceModel( const ArrayType& delta, const ArrayType& K, + const ArrayType& _G0, const ParticleType& _type ) + : base_type( delta, K, _type ) + , G0( view_type_1d( "fracture_energy", num_types ) ) + , s0( view_type_2d( "critical_stretch", num_types, num_types ) ) + , bond_break_coeff( + view_type_2d( "break_coeff", num_types, num_types ) ) + { + setParameters( _G0 ); + } + + template + void setParameters( const ArrayType& _G0 ) + { + // Initialize self interaction parameters. + auto init_self_func = KOKKOS_LAMBDA( const int i ) + { + G0( i ) = _G0[i]; + s0( i, i ) = criticalStretch( i ); + bond_break_coeff( i, i ) = + ( 1.0 + s0( i, i ) ) * ( 1.0 + s0( i, i ) ); + }; + using exec_space = typename memory_space::execution_space; + Kokkos::RangePolicy policy( 0, num_types ); + Kokkos::parallel_for( "CabanaPD::Model::Init", policy, init_self_func ); + Kokkos::fence(); + + // Initialize cross-terms. + auto init_cross_func = KOKKOS_LAMBDA( const int i ) + { + for ( std::size_t j = i; j < num_types; j++ ) + { + s0( i, j ) = criticalStretch( i, j ); + bond_break_coeff( i, j ) = + ( 1.0 + s0( i, j ) ) * ( 1.0 + s0( i, j ) ); + } + }; + Kokkos::parallel_for( "CabanaPD::Model::Init", policy, + init_cross_func ); + } + + KOKKOS_INLINE_FUNCTION + auto criticalStretch( const int type_i ) const + { + return sqrt( 5.0 * G0( type_i ) / 9.0 / K( type_i ) / delta( type_i ) ); + } + + KOKKOS_INLINE_FUNCTION + auto criticalStretch( const int type_i, const int type_j ) const + { + auto s0_i = s0( type_i, type_i ); + auto s0_j = s0( type_j, type_j ); + auto c_i = c( type_i, type_i ); + auto c_j = c( type_j, type_j ); + return Kokkos::sqrt( ( s0_i * s0_i * c_i + s0_j * s0_j * c_j ) / + ( c_i + c_j ) ); + } + + KOKKOS_INLINE_FUNCTION + bool criticalStretch( const int i, const int j, const double r, + const double xi ) const + { + // if ( type( j ) != 0 ) + // std::cout << bond_break_coeff.size() << " " << type( i ) << " " + // << type( j ) << "\n"; + //<< type.size() << " " << i << " " << j << "\n"; + return r * r >= bond_break_coeff( type( i ), type( j ) ) * xi * xi; + } +}; + +template +auto createForceModel( PMB, Fracture, TemperatureIndependent, + ParticleType particles, const ArrayType& delta, + const ArrayType& K, const ArrayType& G0 ) +{ + auto type = particles.sliceType(); + using type_type = decltype( type ); + return ForceModel( + delta, K, G0, type ); +} + +template +struct ForceModel + : public ForceModel +{ + using base_type = + ForceModel; using base_model = typename base_type::base_model; + using species_type = typename base_type::species_type; using fracture_type = typename base_type::fracture_type; - using thermal_type = base_type::thermal_type; + using thermal_type = typename base_type::thermal_type; using base_type::base_type; @@ -87,14 +302,16 @@ struct ForceModel using base_type::K; }; -template <> -struct ForceModel - : public ForceModel +template +struct ForceModel + : public ForceModel { - using base_type = ForceModel; + using base_type = + ForceModel; using base_model = typename base_type::base_model; + using species_type = typename base_type::species_type; using fracture_type = typename base_type::fracture_type; - using thermal_type = base_type::thermal_type; + using thermal_type = typename base_type::thermal_type; using base_type::base_type; @@ -112,7 +329,9 @@ struct ForceModel : public ForceModel, BaseTemperatureModel { + using memory_space = typename TemperatureType::memory_space; using base_type = ForceModel; + using species_type = typename base_type::species_type; using base_temperature_type = BaseTemperatureModel; using base_model = PMB; using fracture_type = Elastic; @@ -140,7 +359,8 @@ struct ForceModel }; template -auto createForceModel( PMB, Elastic, ParticleType particles, const double delta, +auto createForceModel( PMB, Elastic, TemperatureDependent, + ParticleType particles, const double delta, const double K, const double alpha, const double temp0 ) { auto temp = particles.sliceTemperature(); @@ -154,7 +374,9 @@ struct ForceModel : public ForceModel, BaseTemperatureModel { + using memory_space = typename TemperatureType::memory_space; using base_type = ForceModel; + using species_type = typename base_type::species_type; using base_temperature_type = BaseTemperatureModel; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; diff --git a/src/force/CabanaPD_Force_LPS.hpp b/src/force/CabanaPD_Force_LPS.hpp index 78f25f81..c14d3c2c 100644 --- a/src/force/CabanaPD_Force_LPS.hpp +++ b/src/force/CabanaPD_Force_LPS.hpp @@ -69,21 +69,25 @@ namespace CabanaPD { -template -class Force> +template +class Force> + : public ForceBase { + public: + using exec_space = ExecutionSpace; + using model_type = ForceModel; + using base_type = ForceBase; + protected: - bool _half_neigh; - ForceModel _model; + using base_type::_half_neigh; + model_type _model; - Timer _timer; - Timer _energy_timer; + using base_type::_energy_timer; + using base_type::_timer; public: - using exec_space = ExecutionSpace; - Force( const bool half_neigh, const ForceModel model ) - : _half_neigh( half_neigh ) + : base_type( half_neigh ) , _model( model ) { } @@ -247,9 +251,6 @@ class Force> _energy_timer.stop(); return strain_energy; } - - auto time() { return _timer.time(); }; - auto timeEnergy() { return _energy_timer.time(); }; }; template @@ -498,22 +499,25 @@ class Force> } }; -template -class Force> - : public Force> +template +class Force> + : public Force> { + public: + using exec_space = ExecutionSpace; + using model_type = ForceModel; + using base_type = + Force>; + protected: - using base_type = Force>; using base_type::_half_neigh; - ForceModel _model; + model_type _model; using base_type::_energy_timer; using base_type::_timer; public: - using exec_space = ExecutionSpace; - - Force( const bool half_neigh, const ForceModel model ) + Force( const bool half_neigh, const model_type model ) : base_type( half_neigh, model ) , _model( model ) { diff --git a/src/force/CabanaPD_Force_PMB.hpp b/src/force/CabanaPD_Force_PMB.hpp index 5b8cd537..8c3ecd15 100644 --- a/src/force/CabanaPD_Force_PMB.hpp +++ b/src/force/CabanaPD_Force_PMB.hpp @@ -72,12 +72,12 @@ namespace CabanaPD { template class Force> - : public Force + : public ForceBase { public: using exec_space = ExecutionSpace; using model_type = ForceModel; - using base_type = Force; + using base_type = ForceBase; protected: using base_type::_half_neigh; @@ -117,7 +117,7 @@ class Force> model.thermalStretch( s, i, j ); - const double coeff = model.c * s * vol( j ); + const double coeff = model.micromodulus( i, j ) * s * vol( j ); fx_i = coeff * rx / r; fy_i = coeff * ry / r; fz_i = coeff * rz / r; @@ -158,7 +158,8 @@ class Force> // 0.25 factor is due to 1/2 from outside the integral and 1/2 from // the integrand (pairwise potential). - double w = 0.25 * model.c * s * s * xi * vol( j ); + double w = + 0.25 * model.micromodulus( i, j ) * s * s * xi * vol( j ); W( i ) += w; Phi += w * vol( i ); }; @@ -177,7 +178,7 @@ class Force> template class Force> - : public Force + : public ForceBase { public: using exec_space = ExecutionSpace; @@ -185,7 +186,7 @@ class Force> protected: using base_model_type = typename model_type::base_type; - using base_type = Force; + using base_type = ForceBase; using base_type::_half_neigh; model_type _model; @@ -243,7 +244,8 @@ class Force> // Else if statement is only for performance. else if ( mu( i, n ) > 0 ) { - const double coeff = model.c * s * vol( j ); + const double coeff = + model.micromodulus( i, j ) * s * vol( j ); double muij = mu( i, n ); fx_i = muij * coeff * rx / r; fy_i = muij * coeff * ry / r; @@ -295,7 +297,8 @@ class Force> // 0.25 factor is due to 1/2 from outside the integral and 1/2 // from the integrand (pairwise potential). - double w = mu( i, n ) * 0.25 * model.c * s * s * xi * vol( j ); + double w = mu( i, n ) * 0.25 * model.micromodulus( i, j ) * s * + s * xi * vol( j ); W( i ) += w; phi_i += mu( i, n ) * vol( j ); @@ -317,14 +320,14 @@ class Force> template class Force> - : public Force + : public ForceBase { public: using exec_space = ExecutionSpace; - using model_type = ForceModel; + using model_type = ForceModel; protected: - using base_type = Force; + using base_type = ForceBase; using base_type::_half_neigh; model_type _model; diff --git a/unit_test/tstComm.hpp b/unit_test/tstComm.hpp index bcbbceea..de45a8b6 100644 --- a/unit_test/tstComm.hpp +++ b/unit_test/tstComm.hpp @@ -77,7 +77,7 @@ void testHalo() Cabana::deep_copy( rank_init_host, rank ); // A gather is performed on construction. - CabanaPD::Comm comm( particles ); From dc0101c3425da990c327ba4809ba313c0532530b Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 17 Jul 2024 15:38:22 -0400 Subject: [PATCH 3/7] DO NOT MERGE; test with KW --- examples/mechanics/kalthoff_winkler.cpp | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/examples/mechanics/kalthoff_winkler.cpp b/examples/mechanics/kalthoff_winkler.cpp index 4ed9159c..08e8c120 100644 --- a/examples/mechanics/kalthoff_winkler.cpp +++ b/examples/mechanics/kalthoff_winkler.cpp @@ -76,23 +76,26 @@ void kalthoffWinklerExample( const std::string filename ) CabanaPD::Prenotch<2> prenotch( v1, v2, notch_positions ); // ==================================================== - // Force model + // Force model options // ==================================================== - using model_type = CabanaPD::ForceModel; - model_type force_model( delta, K, G0 ); - // using model_type = - // CabanaPD::ForceModel; - // model_type force_model( delta, K, G, G0 ); + using model_type = CabanaPD::PMB; + using thermal_type = CabanaPD::TemperatureIndependent; // ==================================================== // Particle generation // ==================================================== // Does not set displacements, velocities, etc. auto particles = std::make_shared< - CabanaPD::Particles>( + CabanaPD::Particles>( exec_space(), low_corner, high_corner, num_cells, halo_width ); + std::array same_d = { delta, delta }; + std::array fake_K = { K, K * 1.2 }; + std::array fake_G0 = { G0, G0 * 1.8 }; + auto force_model = + createForceModel( model_type{}, CabanaPD::Fracture{}, thermal_type{}, + *particles, same_d, fake_K, fake_G0 ); + // ==================================================== // Boundary conditions // ==================================================== @@ -111,6 +114,7 @@ void kalthoffWinklerExample( const std::string filename ) auto x = particles->sliceReferencePosition(); auto v = particles->sliceVelocity(); auto f = particles->sliceForce(); + auto type = particles->sliceType(); double v0 = inputs["impactor_velocity"]; auto init_functor = KOKKOS_LAMBDA( const int pid ) @@ -121,6 +125,8 @@ void kalthoffWinklerExample( const std::string filename ) if ( x( pid, 1 ) > y_prenotch1 && x( pid, 1 ) < y_prenotch2 && x( pid, 0 ) < -0.5 * height + dx ) v( pid, 0 ) = v0; + if ( x( pid, 1 ) > 0.0 ) + type( pid ) = 1; }; particles->updateParticles( exec_space{}, init_functor ); From a8ce5c4dfaa427f6dc12dbd113793bee096c80df Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Tue, 23 Jul 2024 12:30:31 -0400 Subject: [PATCH 4/7] Require c++ 20 and use class lambda for force models --- .github/workflows/CI.yml | 3 +++ src/CabanaPD_ForceModels.hpp | 2 +- src/force/CabanaPD_ForceModels_PMB.hpp | 8 ++++---- 3 files changed, 8 insertions(+), 5 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 01f68736..4420f25f 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -57,6 +57,7 @@ jobs: -DCMAKE_BUILD_TYPE=${{ matrix.cmake_build_type }} \ -DCMAKE_CXX_COMPILER=${{ matrix.cxx }} \ -DCMAKE_INSTALL_PREFIX=$HOME/kokkos \ + -DCMAKE_CXX_STANDARD=20 \ -DKokkos_ENABLE_${{ matrix.backend }}=ON \ -DKokkos_ENABLE_HWLOC=ON cmake --build build --parallel 2 @@ -148,6 +149,7 @@ jobs: run: | cmake -B build \ -DCMAKE_INSTALL_PREFIX=$HOME/kokkos \ + -DCMAKE_CXX_STANDARD=20 \ -DKokkos_ENABLE_HIP=ON \ -DKokkos_ARCH_VEGA908=ON \ -DCMAKE_BUILD_TYPE=${{ matrix.cmake_build_type }} \ @@ -230,6 +232,7 @@ jobs: run: | cmake -B build \ -DCMAKE_INSTALL_PREFIX=$HOME/kokkos \ + -DCMAKE_CXX_STANDARD=20 \ -DKokkos_ENABLE_CUDA=ON \ -DKokkos_ARCH_VOLTA72=ON \ -DKokkos_ENABLE_CUDA_LAMBDA=ON \ diff --git a/src/CabanaPD_ForceModels.hpp b/src/CabanaPD_ForceModels.hpp index 93afb9ca..a4b9c6d5 100644 --- a/src/CabanaPD_ForceModels.hpp +++ b/src/CabanaPD_ForceModels.hpp @@ -66,7 +66,7 @@ struct BaseForceModel , num_types( _delta.size() ) { max_delta = 0; - auto init_func = KOKKOS_LAMBDA( const int i, double& max ) + auto init_func = KOKKOS_CLASS_LAMBDA( const int i, double& max ) { delta( i ) = _delta[i]; if ( delta( i ) > max ) diff --git a/src/force/CabanaPD_ForceModels_PMB.hpp b/src/force/CabanaPD_ForceModels_PMB.hpp index 915fa9d0..dbb6fe0c 100644 --- a/src/force/CabanaPD_ForceModels_PMB.hpp +++ b/src/force/CabanaPD_ForceModels_PMB.hpp @@ -92,7 +92,7 @@ struct ForceModel void setParameters( const ArrayType& _K ) { // Initialize self interaction parameters. - auto init_self_func = KOKKOS_LAMBDA( const int i ) + auto init_self_func = KOKKOS_CLASS_LAMBDA( const int i ) { K( i ) = _K[i]; c( i, i ) = micromodulus( i ); @@ -103,7 +103,7 @@ struct ForceModel Kokkos::fence(); // Initialize cross-terms. - auto init_cross_func = KOKKOS_LAMBDA( const int i ) + auto init_cross_func = KOKKOS_CLASS_LAMBDA( const int i ) { for ( std::size_t j = i; j < num_types; j++ ) c( i, j ) = ( micromodulus( i ) + micromodulus( j ) ) / 2.0; @@ -218,7 +218,7 @@ struct ForceModel void setParameters( const ArrayType& _G0 ) { // Initialize self interaction parameters. - auto init_self_func = KOKKOS_LAMBDA( const int i ) + auto init_self_func = KOKKOS_CLASS_LAMBDA( const int i ) { G0( i ) = _G0[i]; s0( i, i ) = criticalStretch( i ); @@ -231,7 +231,7 @@ struct ForceModel Kokkos::fence(); // Initialize cross-terms. - auto init_cross_func = KOKKOS_LAMBDA( const int i ) + auto init_cross_func = KOKKOS_CLASS_LAMBDA( const int i ) { for ( std::size_t j = i; j < num_types; j++ ) { From 8ff96ed477aa14ff5fcf71694983fca747027e16 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 9 Aug 2024 11:34:09 -0400 Subject: [PATCH 5/7] Add multi-species thermal models; simplifies template parameters with single-species tag --- src/CabanaPD_ForceModels.hpp | 61 ++++++++++++-- src/force/CabanaPD_ForceModels_LPS.hpp | 4 +- src/force/CabanaPD_ForceModels_PMB.hpp | 110 +++++++++++++++++++++++-- 3 files changed, 158 insertions(+), 17 deletions(-) diff --git a/src/CabanaPD_ForceModels.hpp b/src/CabanaPD_ForceModels.hpp index a4b9c6d5..a75b7e41 100644 --- a/src/CabanaPD_ForceModels.hpp +++ b/src/CabanaPD_ForceModels.hpp @@ -17,11 +17,11 @@ namespace CabanaPD { -template +template struct BaseForceModel; template <> -struct BaseForceModel<> +struct BaseForceModel { using species_type = SingleSpecies; double delta; @@ -38,14 +38,13 @@ struct BaseForceModel<> void thermalStretch( double&, const int, const int ) const {} }; -template +template struct BaseForceModel { using species_type = MultiSpecies; // Only allow one memory space. - using memory_space = - typename std::tuple_element<0, std::tuple>::type; + using memory_space = MemorySpace; using view_type_1d = Kokkos::View; view_type_1d delta; double max_delta; @@ -86,8 +85,11 @@ struct BaseForceModel void thermalStretch( double&, const int, const int ) const {} }; +template +struct BaseTemperatureModel; + template -struct BaseTemperatureModel +struct BaseTemperatureModel { using species_type = SingleSpecies; using memory_space = typename TemperatureType::memory_space; @@ -115,6 +117,53 @@ struct BaseTemperatureModel } }; +template +struct BaseTemperatureModel +{ + using species_type = MultiSpecies; + using memory_space = typename TemperatureType::memory_space; + using view_type_1d = Kokkos::View; + view_type_1d alpha; + view_type_1d temp0; + + // Particle fields + TemperatureType temperature; + ParticleType type; + + template + BaseTemperatureModel( const TemperatureType& _temp, const ArrayType _alpha, + const ArrayType _temp0, const ParticleType& _type ) + : alpha( view_type_1d( "delta", _alpha.size() ) ) + , temp0( view_type_1d( "delta", _temp0.size() ) ) + , temperature( _temp ) + , type( _type ) + { + auto init_func = KOKKOS_CLASS_LAMBDA( const int i ) + { + alpha( i ) = _alpha[i]; + temp0( i ) = _temp0[i]; + }; + using exec_space = typename memory_space::execution_space; + Kokkos::RangePolicy policy( 0, alpha.size() ); + Kokkos::parallel_for( "CabanaPD::Model::Init", policy, init_func ); + } + + void update( const TemperatureType& _temp, const ParticleType& _type ) + { + temperature = _temp; + type = _type; + } + + // Update stretch with temperature effects. + KOKKOS_INLINE_FUNCTION + void thermalStretch( double& s, const int i, const int j ) const + { + double temp_avg = + 0.5 * ( temperature( i ) + temperature( j ) ) - temp0( type( i ) ); + s -= alpha( type( i ) ) * temp_avg; + } +}; + template struct ForceModel; diff --git a/src/force/CabanaPD_ForceModels_LPS.hpp b/src/force/CabanaPD_ForceModels_LPS.hpp index 14bd9be8..b42eb30c 100644 --- a/src/force/CabanaPD_ForceModels_LPS.hpp +++ b/src/force/CabanaPD_ForceModels_LPS.hpp @@ -19,9 +19,9 @@ namespace CabanaPD { template <> -struct ForceModel : public BaseForceModel<> +struct ForceModel : public BaseForceModel { - using base_type = BaseForceModel<>; + using base_type = BaseForceModel; using species_type = SingleSpecies; using base_model = LPS; using fracture_type = Elastic; diff --git a/src/force/CabanaPD_ForceModels_PMB.hpp b/src/force/CabanaPD_ForceModels_PMB.hpp index dbb6fe0c..55d189ce 100644 --- a/src/force/CabanaPD_ForceModels_PMB.hpp +++ b/src/force/CabanaPD_ForceModels_PMB.hpp @@ -21,10 +21,10 @@ namespace CabanaPD template <> struct ForceModel - : public BaseForceModel<> + : public BaseForceModel { - using base_type = BaseForceModel<>; - using species_type = typename base_type::species_type; + using base_type = BaseForceModel; + using species_type = SingleSpecies; using base_model = PMB; using fracture_type = Elastic; using thermal_type = TemperatureIndependent; @@ -327,12 +327,13 @@ struct ForceModel template struct ForceModel : public ForceModel, - BaseTemperatureModel + BaseTemperatureModel { using memory_space = typename TemperatureType::memory_space; using base_type = ForceModel; using species_type = typename base_type::species_type; - using base_temperature_type = BaseTemperatureModel; + using base_temperature_type = + BaseTemperatureModel; using base_model = PMB; using fracture_type = Elastic; using thermal_type = TemperatureDependent; @@ -348,7 +349,6 @@ struct ForceModel // Explicitly use the temperature-dependent stretch. using base_temperature_type::thermalStretch; - // ForceModel(){}; ForceModel( const double _delta, const double _K, const TemperatureType _temp, const double _alpha, const double _temp0 = 0.0 ) @@ -358,6 +358,44 @@ struct ForceModel } }; +template +struct ForceModel + : public ForceModel, + BaseTemperatureModel +{ + using memory_space = typename TemperatureType::memory_space; + using base_type = + ForceModel; + using species_type = typename base_type::species_type; + using base_temperature_type = + BaseTemperatureModel; + using base_model = PMB; + using fracture_type = Elastic; + using thermal_type = TemperatureDependent; + + using base_type::c; + using base_type::delta; + using base_type::K; + using base_type::type; + + // Thermal parameters + using base_temperature_type::alpha; + using base_temperature_type::temp0; + + // Explicitly use the temperature-dependent stretch. + using base_temperature_type::thermalStretch; + + template + ForceModel( const ArrayType _delta, const ArrayType _K, + const TemperatureType& _temp, const ArrayType _alpha, + const ArrayType _temp0, ParticleType& _type ) + : base_type( _delta, _K, _type ) + , base_temperature_type( _temp, _alpha, _temp0 ) + { + } +}; + template auto createForceModel( PMB, Elastic, TemperatureDependent, ParticleType particles, const double delta, @@ -372,12 +410,13 @@ auto createForceModel( PMB, Elastic, TemperatureDependent, template struct ForceModel : public ForceModel, - BaseTemperatureModel + BaseTemperatureModel { using memory_space = typename TemperatureType::memory_space; using base_type = ForceModel; using species_type = typename base_type::species_type; - using base_temperature_type = BaseTemperatureModel; + using base_temperature_type = + BaseTemperatureModel; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; using thermal_type = TemperatureDependent; @@ -398,7 +437,6 @@ struct ForceModel // Explicitly use the temperature-dependent stretch. using base_temperature_type::thermalStretch; - // ForceModel(){}; ForceModel( const double _delta, const double _K, const double _G0, const TemperatureType _temp, const double _alpha, const double _temp0 = 0.0 ) @@ -418,6 +456,60 @@ struct ForceModel } }; +template +struct ForceModel + : public ForceModel, + BaseTemperatureModel +{ + using memory_space = typename TemperatureType::memory_space; + using base_type = + ForceModel; + using species_type = typename base_type::species_type; + using base_temperature_type = + BaseTemperatureModel; + using base_model = typename base_type::base_model; + using fracture_type = typename base_type::fracture_type; + using thermal_type = TemperatureDependent; + + using base_type::c; + using base_type::delta; + using base_type::K; + using base_type::type; + + // Does not use the base bond_break_coeff. + using base_type::G0; + using base_type::s0; + + // Thermal parameters + using base_temperature_type::alpha; + using base_temperature_type::temp0; + using base_temperature_type::temperature; + + // Explicitly use the temperature-dependent stretch. + using base_temperature_type::thermalStretch; + + template + ForceModel( const ArrayType _delta, const ArrayType _K, const ArrayType _G0, + const TemperatureType _temp, const ArrayType _alpha, + const ArrayType _temp0, ParticleType& _type ) + : base_type( _delta, _K, _G0, _type ) + , base_temperature_type( _temp, _alpha, _temp0 ) + { + } + + KOKKOS_INLINE_FUNCTION + bool criticalStretch( const int i, const int j, const double r, + const double xi ) const + { + double temp_avg = + 0.5 * ( temperature( i ) + temperature( j ) ) - temp0( type( i ) ); + double bond_break_coeff = ( 1.0 + s0 + alpha( type( i ) ) * temp_avg ) * + ( 1.0 + s0 + alpha( type( i ) ) * temp_avg ); + return r * r >= bond_break_coeff * xi * xi; + } +}; + template auto createForceModel( PMB, Fracture, ParticleType particles, const double delta, const double K, const double G0, From c089abf5445f5124635d0293c0dd6c7e16ea3272 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Tue, 20 Aug 2024 12:58:49 -0400 Subject: [PATCH 6/7] fixup: remove class lambdas; use local copies + class operators --- .github/workflows/CI.yml | 3 -- src/CabanaPD_ForceModels.hpp | 34 ++++++++++---- src/force/CabanaPD_ForceModels_PMB.hpp | 65 +++++++++++++++----------- 3 files changed, 61 insertions(+), 41 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 4420f25f..01f68736 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -57,7 +57,6 @@ jobs: -DCMAKE_BUILD_TYPE=${{ matrix.cmake_build_type }} \ -DCMAKE_CXX_COMPILER=${{ matrix.cxx }} \ -DCMAKE_INSTALL_PREFIX=$HOME/kokkos \ - -DCMAKE_CXX_STANDARD=20 \ -DKokkos_ENABLE_${{ matrix.backend }}=ON \ -DKokkos_ENABLE_HWLOC=ON cmake --build build --parallel 2 @@ -149,7 +148,6 @@ jobs: run: | cmake -B build \ -DCMAKE_INSTALL_PREFIX=$HOME/kokkos \ - -DCMAKE_CXX_STANDARD=20 \ -DKokkos_ENABLE_HIP=ON \ -DKokkos_ARCH_VEGA908=ON \ -DCMAKE_BUILD_TYPE=${{ matrix.cmake_build_type }} \ @@ -232,7 +230,6 @@ jobs: run: | cmake -B build \ -DCMAKE_INSTALL_PREFIX=$HOME/kokkos \ - -DCMAKE_CXX_STANDARD=20 \ -DKokkos_ENABLE_CUDA=ON \ -DKokkos_ARCH_VOLTA72=ON \ -DKokkos_ENABLE_CUDA_LAMBDA=ON \ diff --git a/src/CabanaPD_ForceModels.hpp b/src/CabanaPD_ForceModels.hpp index a75b7e41..9a1ff501 100644 --- a/src/CabanaPD_ForceModels.hpp +++ b/src/CabanaPD_ForceModels.hpp @@ -63,13 +63,20 @@ struct BaseForceModel BaseForceModel( const ArrayType& _delta ) : delta( view_type_1d( "delta", _delta.size() ) ) , num_types( _delta.size() ) + { + setParameters( _delta ); + } + + template + void setParameters( const ArrayType& _delta ) { max_delta = 0; - auto init_func = KOKKOS_CLASS_LAMBDA( const int i, double& max ) + auto delta_copy = delta; + auto init_func = KOKKOS_LAMBDA( const int i, double& max ) { - delta( i ) = _delta[i]; - if ( delta( i ) > max ) - max = delta( i ); + delta_copy( i ) = _delta[i]; + if ( delta_copy( i ) > max ) + max = delta_copy( i ); }; using exec_space = typename memory_space::execution_space; Kokkos::RangePolicy policy( 0, num_types ); @@ -99,7 +106,6 @@ struct BaseTemperatureModel // Temperature field TemperatureType temperature; - BaseTemperatureModel(){}; BaseTemperatureModel( const TemperatureType _temp, const double _alpha, const double _temp0 ) : alpha( _alpha ) @@ -131,17 +137,25 @@ struct BaseTemperatureModel ParticleType type; template - BaseTemperatureModel( const TemperatureType& _temp, const ArrayType _alpha, - const ArrayType _temp0, const ParticleType& _type ) + BaseTemperatureModel( const TemperatureType& _temp, const ArrayType& _alpha, + const ArrayType& _temp0, const ParticleType& _type ) : alpha( view_type_1d( "delta", _alpha.size() ) ) , temp0( view_type_1d( "delta", _temp0.size() ) ) , temperature( _temp ) , type( _type ) { - auto init_func = KOKKOS_CLASS_LAMBDA( const int i ) + setParameters( _alpha, _temp0 ); + } + + template + void setParameters( const ArrayType& _alpha, const ArrayType& _temp0 ) + { + auto alpha_copy = alpha; + auto temp0_copy = temp0; + auto init_func = KOKKOS_LAMBDA( const int i ) { - alpha( i ) = _alpha[i]; - temp0( i ) = _temp0[i]; + alpha_copy( i ) = _alpha[i]; + temp0_copy( i ) = _temp0[i]; }; using exec_space = typename memory_space::execution_space; Kokkos::RangePolicy policy( 0, alpha.size() ); diff --git a/src/force/CabanaPD_ForceModels_PMB.hpp b/src/force/CabanaPD_ForceModels_PMB.hpp index 55d189ce..dd256aeb 100644 --- a/src/force/CabanaPD_ForceModels_PMB.hpp +++ b/src/force/CabanaPD_ForceModels_PMB.hpp @@ -91,25 +91,30 @@ struct ForceModel template void setParameters( const ArrayType& _K ) { - // Initialize self interaction parameters. - auto init_self_func = KOKKOS_CLASS_LAMBDA( const int i ) + // Initialize per-type variables. + auto K_copy = K; + auto init_self_func = KOKKOS_LAMBDA( const int i ) { - K( i ) = _K[i]; - c( i, i ) = micromodulus( i ); + K_copy( i ) = _K[i]; }; using exec_space = typename memory_space::execution_space; Kokkos::RangePolicy policy( 0, num_types ); - Kokkos::parallel_for( "CabanaPD::Model::Init", policy, init_self_func ); + Kokkos::parallel_for( "CabanaPD::Model::Copy", policy, init_self_func ); Kokkos::fence(); - // Initialize cross-terms. - auto init_cross_func = KOKKOS_CLASS_LAMBDA( const int i ) + // Initialize model parameters. + Kokkos::parallel_for( "CabanaPD::Model::Init", policy, *this ); + } + + KOKKOS_INLINE_FUNCTION void operator()( const int i ) const + { + c( i, i ) = micromodulus( i ); + for ( std::size_t j = i; j < num_types; j++ ) { - for ( std::size_t j = i; j < num_types; j++ ) - c( i, j ) = ( micromodulus( i ) + micromodulus( j ) ) / 2.0; - }; - Kokkos::parallel_for( "CabanaPD::Model::Init", policy, - init_cross_func ); + c( i, j ) = ( micromodulus( i ) + micromodulus( j ) ) / 2.0; + // Set symmetric cross-terms. + c( j, i ) = c( i, j ); + } } KOKKOS_INLINE_FUNCTION @@ -217,13 +222,11 @@ struct ForceModel template void setParameters( const ArrayType& _G0 ) { - // Initialize self interaction parameters. - auto init_self_func = KOKKOS_CLASS_LAMBDA( const int i ) + // Initialize per-type variables. + auto G0_copy = G0; + auto init_self_func = KOKKOS_LAMBDA( const int i ) { - G0( i ) = _G0[i]; - s0( i, i ) = criticalStretch( i ); - bond_break_coeff( i, i ) = - ( 1.0 + s0( i, i ) ) * ( 1.0 + s0( i, i ) ); + G0_copy( i ) = _G0[i]; }; using exec_space = typename memory_space::execution_space; Kokkos::RangePolicy policy( 0, num_types ); @@ -231,17 +234,23 @@ struct ForceModel Kokkos::fence(); // Initialize cross-terms. - auto init_cross_func = KOKKOS_CLASS_LAMBDA( const int i ) + Kokkos::parallel_for( "CabanaPD::Model::Init", policy, *this ); + } + + KOKKOS_INLINE_FUNCTION void operator()( const int i ) const + { + s0( i, i ) = criticalStretch( i ); + bond_break_coeff( i, i ) = ( 1.0 + s0( i, i ) ) * ( 1.0 + s0( i, i ) ); + + for ( std::size_t j = i; j < num_types; j++ ) { - for ( std::size_t j = i; j < num_types; j++ ) - { - s0( i, j ) = criticalStretch( i, j ); - bond_break_coeff( i, j ) = - ( 1.0 + s0( i, j ) ) * ( 1.0 + s0( i, j ) ); - } - }; - Kokkos::parallel_for( "CabanaPD::Model::Init", policy, - init_cross_func ); + s0( i, j ) = criticalStretch( i, j ); + bond_break_coeff( i, j ) = + ( 1.0 + s0( i, j ) ) * ( 1.0 + s0( i, j ) ); + // Set symmetric cross-terms. + s0( j, i ) = s0( i, j ); + bond_break_coeff( j, i ) = bond_break_coeff( i, j ); + } } KOKKOS_INLINE_FUNCTION From 041f7bca6976c542c36cde2c774e2c1e3e3a5546 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 30 Oct 2024 14:46:49 -0400 Subject: [PATCH 7/7] Very WIP --- .../inputs/thermal_deformation.json | 3 +- src/CabanaPD_ForceModels.hpp | 1 - src/CabanaPD_Properties.hpp | 109 ++++++++++++++++++ src/force/CabanaPD_ForceModels_PMB.hpp | 30 +++-- 4 files changed, 134 insertions(+), 9 deletions(-) create mode 100644 src/CabanaPD_Properties.hpp diff --git a/examples/thermomechanics/inputs/thermal_deformation.json b/examples/thermomechanics/inputs/thermal_deformation.json index 4f0f6908..fcca5233 100644 --- a/examples/thermomechanics/inputs/thermal_deformation.json +++ b/examples/thermomechanics/inputs/thermal_deformation.json @@ -1,8 +1,9 @@ { "num_cells" : {"value": [100, 30, 3]}, "system_size" : {"value": [1.0, 0.3, 0.03], "unit": "m"}, + "num_types" : {"value": 1}, "density" : {"value": 3980, "unit": "kg/m^3"}, - "elastic_modulus" : {"value": 370e+9, "unit": "Pa"}, + "elastic_modulus" : {"value": [370e+9, 100e+9, 10e+9], "unit": "Pa", "note": "11, 22, 12"}, "thermal_coefficient" : {"value": 7.5E-6, "unit": "oC^{-1}"}, "reference_temperature" : {"value": 20.0, "unit": "oC"}, "horizon" : {"value": 0.03, "unit": "m"}, diff --git a/src/CabanaPD_ForceModels.hpp b/src/CabanaPD_ForceModels.hpp index 9a1ff501..e5f6cd86 100644 --- a/src/CabanaPD_ForceModels.hpp +++ b/src/CabanaPD_ForceModels.hpp @@ -43,7 +43,6 @@ struct BaseForceModel { using species_type = MultiSpecies; - // Only allow one memory space. using memory_space = MemorySpace; using view_type_1d = Kokkos::View; view_type_1d delta; diff --git a/src/CabanaPD_Properties.hpp b/src/CabanaPD_Properties.hpp new file mode 100644 index 00000000..c5b2979b --- /dev/null +++ b/src/CabanaPD_Properties.hpp @@ -0,0 +1,109 @@ +/**************************************************************************** + * 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 PROPERTIES_HPP +#define PROPERTIES_HPP + +#include + +namespace CabanaPD +{ +//---------------------------------------------------------------------------// +// Properties. +//---------------------------------------------------------------------------// +namespace Field +{ + +template +class Property; + +template <> +class Property +{ + using species_type = SingleSpecies; + double value; + + Property( const double _value ) + : value( _value ){}; +}; + +template +struct Property +{ + using species_type = MultiSpecies; + + using memory_space = MemorySpace; + using view_type_1d = Kokkos::View; + view_type_1d delta; + double max_delta; + std::size_t num_types; + + using view_type_2d = Kokkos::View; + + using view_type_map = Kokkos::View; + view_type_map map; + + Property( const double _delta ) + : delta( view_type_1d( "delta", 1 ) ) + , num_types( 1 ) + { + Kokkos::deep_copy( delta, _delta ); + } + + template + Property( const int _num_types, const ArrayType& _delta ) + : delta( view_type_1d( "delta", _delta.size() ) ) + , num_types( _num_types ) + , map( view_type_map( "csr_to_2d_indexing", _delta.size() ) ) + { + // 00 01 02 03 0 5 4 0 2 + // 10 11 12 13 1 3 1 + // 20 21 22 23 2 + // 30 31 32 33 + std::vector vec_i; + std::vector vec_j; + for ( std::size_t j = i; j < num_types; j++ ) + vec_i[i] = ; + // i 0 1 2 3 2 1 0 0 0 1 + // j 0 1 2 3 3 3 3 2 1 2 + // 0 1 2 1 0 0 + // 0 1 0 + setParameters( _delta ); + } + + template + void setParameters( const ArrayType& _K ) + { + // Copy into Views. + auto K_copy = K; + auto num_types_copy = num_types; + auto init_self_func = KOKKOS_LAMBDA( const int i ) + { + K_copy( map_i( i ), map_j( i ) ) = _K[i]; + }; + + using exec_space = typename memory_space::execution_space; + // This could be number of types (averaging for cross terms) + // or number of types and cross terms + Kokkos::RangePolicy policy( 0, _K.size() ); + Kokkos::parallel_for( "CabanaPD::Model::Copy", policy, init_self_func ); + Kokkos::fence(); + + // Initialize model parameters. + if ( _K.size() == num_types ) + Kokkos::parallel_for( "CabanaPD::Model::Init", policy, *this ); + } +}; + +} // namespace Field +} // namespace CabanaPD + +#endif diff --git a/src/force/CabanaPD_ForceModels_PMB.hpp b/src/force/CabanaPD_ForceModels_PMB.hpp index dd256aeb..8fd32c47 100644 --- a/src/force/CabanaPD_ForceModels_PMB.hpp +++ b/src/force/CabanaPD_ForceModels_PMB.hpp @@ -71,10 +71,9 @@ struct ForceModel using base_type::delta; using base_type::num_types; - using view_type_1d = typename base_type::view_type_1d; using view_type_2d = typename base_type::view_type_2d; view_type_2d c; - view_type_1d K; + view_type_2d K; ParticleType type; template @@ -82,7 +81,7 @@ struct ForceModel ParticleType _type ) : base_type( delta ) , c( view_type_2d( "micromodulus", num_types, num_types ) ) - , K( view_type_1d( "bulk_modulus", num_types ) ) + , K( view_type_2d( "bulk_modulus", num_types, num_types ) ) , type( _type ) { setParameters( _K ); @@ -93,22 +92,39 @@ struct ForceModel { // Initialize per-type variables. auto K_copy = K; + auto num_types_copy = num_types; auto init_self_func = KOKKOS_LAMBDA( const int i ) { - K_copy( i ) = _K[i]; + // Diagonal first. + if ( i < num_types_copy ) + { + K_copy( i, i ) = _K[i]; + } + else + { + } }; + // 00 01 02 0 5 4 0 2 + // 10 11 12 1 3 1 + // 20 21 22 2 using exec_space = typename memory_space::execution_space; - Kokkos::RangePolicy policy( 0, num_types ); + // This could be number of types (averaging for cross terms) + // or number of types and cross terms + Kokkos::RangePolicy policy( 0, _K.size() ); Kokkos::parallel_for( "CabanaPD::Model::Copy", policy, init_self_func ); Kokkos::fence(); // Initialize model parameters. - Kokkos::parallel_for( "CabanaPD::Model::Init", policy, *this ); + if ( _K.size() == num_types ) + Kokkos::parallel_for( "CabanaPD::Model::Init", policy, *this ); + else // if (_K.size() == (num_types) * ) + Kokkos::parallel_for( "CabanaPD::Model::Init", policy, *this ); } KOKKOS_INLINE_FUNCTION void operator()( const int i ) const { - c( i, i ) = micromodulus( i ); + if () + c( i, i ) = micromodulus( i ); for ( std::size_t j = i; j < num_types; j++ ) { c( i, j ) = ( micromodulus( i ) + micromodulus( j ) ) / 2.0;