From b1cae3acdf222f02765020c604c09fcaf9dc53b8 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 6 May 2024 09:21:17 -0400 Subject: [PATCH] Add temperature-based models --- src/CabanaPD_ForceModels.hpp | 45 ++++++++++++++++- src/CabanaPD_Types.hpp | 6 +++ src/force/CabanaPD_ForceModels_PMB.hpp | 62 +++++++++++++++++++++-- src/force/CabanaPD_Force_PMB.hpp | 68 ++++++++++++++++---------- 4 files changed, 149 insertions(+), 32 deletions(-) diff --git a/src/CabanaPD_ForceModels.hpp b/src/CabanaPD_ForceModels.hpp index af1cf283..88bf7be5 100644 --- a/src/CabanaPD_ForceModels.hpp +++ b/src/CabanaPD_ForceModels.hpp @@ -23,9 +23,52 @@ struct BaseForceModel : delta( _delta ){}; BaseForceModel( const BaseForceModel& model ) { delta = model.delta; } + + // No-op for temperature. + KOKKOS_INLINE_FUNCTION + void thermalStretch( double&, const int, const int ) const {} +}; + +template +struct BaseTemperatureModel +{ + double alpha; + double temp0; + + // Temperature field + TemperatureType temperature; + + BaseTemperatureModel(){}; + BaseTemperatureModel( const TemperatureType _temp, const double _alpha, + const double _temp0 ) + : alpha( _alpha ) + , temp0( _temp0 ) + , temperature( _temp ){}; + + BaseTemperatureModel( const BaseTemperatureModel& model ) + { + alpha = model.alpha; + temp0 = model.temp0; + temperature = model.temperature; + } + + void set_param( const double _alpha, const double _temp0 ) + { + alpha = _alpha; + temp0 = _temp0; + } + + // 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; + s -= alpha * temp_avg; + } }; -template +template struct ForceModel; } // namespace CabanaPD diff --git a/src/CabanaPD_Types.hpp b/src/CabanaPD_Types.hpp index a036d867..fd95375a 100644 --- a/src/CabanaPD_Types.hpp +++ b/src/CabanaPD_Types.hpp @@ -20,6 +20,12 @@ struct Elastic struct Fracture { }; +struct TemperatureDependent +{ +}; +struct TemperatureIndependent +{ +}; struct PMB { diff --git a/src/force/CabanaPD_ForceModels_PMB.hpp b/src/force/CabanaPD_ForceModels_PMB.hpp index 3d36740a..6f3cbe97 100644 --- a/src/force/CabanaPD_ForceModels_PMB.hpp +++ b/src/force/CabanaPD_ForceModels_PMB.hpp @@ -18,11 +18,12 @@ namespace CabanaPD { template <> -struct ForceModel : public BaseForceModel +struct ForceModel : public BaseForceModel { using base_type = BaseForceModel; using base_model = PMB; using fracture_type = Elastic; + using thermal_type = TemperatureIndependent; using base_type::delta; @@ -52,7 +53,8 @@ struct ForceModel : public BaseForceModel }; template <> -struct ForceModel : public ForceModel +struct ForceModel + : public ForceModel { using base_type = ForceModel; using base_model = typename base_type::base_model; @@ -90,7 +92,8 @@ struct ForceModel : public ForceModel }; template <> -struct ForceModel : public ForceModel +struct ForceModel + : public ForceModel { using base_type = ForceModel; using base_model = typename base_type::base_model; @@ -104,7 +107,8 @@ struct ForceModel : public ForceModel }; template <> -struct ForceModel : public ForceModel +struct ForceModel + : public ForceModel { using base_type = ForceModel; using base_model = typename base_type::base_model; @@ -121,6 +125,56 @@ struct ForceModel : public ForceModel using base_type::s0; }; +template +struct ForceModel + : public ForceModel, + BaseTemperatureModel +{ + using base_type = ForceModel; + 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; + + // Thermal coefficients + using base_temperature_type::alpha; + using base_temperature_type::temp0; + + // 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 ) + : 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 ); + } +}; + +template +auto createForceModel( ParticleType particles, const double delta, + const double K, const double alpha, const double temp0 ) +{ + auto temp = particles.sliceTemperature(); + using temp_type = decltype( temp ); + return ForceModel( + delta, K, temp, alpha, temp0 ); +} } // namespace CabanaPD #endif diff --git a/src/force/CabanaPD_Force_PMB.hpp b/src/force/CabanaPD_Force_PMB.hpp index 2b1ba5dc..3e53ef19 100644 --- a/src/force/CabanaPD_Force_PMB.hpp +++ b/src/force/CabanaPD_Force_PMB.hpp @@ -69,17 +69,19 @@ namespace CabanaPD { -template -class Force> +template +class Force> { + public: + using exec_space = ExecutionSpace; + using model_type = ForceModel; + protected: bool _half_neigh; - ForceModel _model; + model_type _model; public: - using exec_space = ExecutionSpace; - - Force( const bool half_neigh, const ForceModel model ) + Force( const bool half_neigh, const model_type model ) : _half_neigh( half_neigh ) , _model( model ) { @@ -103,7 +105,7 @@ class Force> const NeighListType& neigh_list, const int n_local, ParallelType& neigh_op_tag ) const { - auto c = _model.c; + auto model = _model; const auto vol = particles.sliceVolume(); auto force_full = KOKKOS_LAMBDA( const int i, const int j ) @@ -115,7 +117,10 @@ class Force> double xi, r, s; double rx, ry, rz; getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz ); - const double coeff = c * s * vol( j ); + + model.thermalStretch( s, i, j ); + + const double coeff = model.c * s * vol( j ); fx_i = coeff * rx / r; fy_i = coeff * ry / r; fz_i = coeff * rz / r; @@ -167,19 +172,22 @@ class Force> } }; -template -class Force> - : public Force> +template +class Force> + : public Force> { + public: + using exec_space = ExecutionSpace; + using model_type = ForceModel; + protected: - using base_type = Force>; + using base_type = + Force>; using base_type::_half_neigh; - ForceModel _model; + model_type _model; 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 ) { @@ -192,7 +200,7 @@ class Force> const NeighListType& neigh_list, MuView& mu, const int n_local, ParallelType& ) const { - auto c = _model.c; + auto model = _model; auto break_coeff = _model.bond_break_coeff; const auto vol = particles.sliceVolume(); const auto nofail = particles.sliceNoFail(); @@ -217,6 +225,8 @@ class Force> double rx, ry, rz; getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz ); + model.thermalStretch( s, i, j ); + // Break if beyond critical stretch unless in no-fail zone. if ( r * r >= break_coeff * xi * xi && !nofail( i ) && !nofail( j ) ) @@ -226,7 +236,7 @@ class Force> // Else if statement is only for performance. else if ( mu( i, n ) > 0 ) { - const double coeff = c * s * vol( j ); + const double coeff = model.c * s * vol( j ); double muij = mu( i, n ); fx_i = muij * coeff * rx / r; fy_i = muij * coeff * ry / r; @@ -291,19 +301,21 @@ class Force> } }; -template -class Force> - : public Force> +template +class Force> + : public Force> { + public: + using exec_space = ExecutionSpace; + using model_type = ForceModel; + protected: using base_type = Force>; using base_type::_half_neigh; - ForceModel _model; + model_type _model; 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 ) { @@ -316,7 +328,7 @@ class Force> const NeighListType& neigh_list, const int n_local, ParallelType& neigh_op_tag ) const { - auto c = _model.c; + auto model = _model; const auto vol = particles.sliceVolume(); auto force_full = KOKKOS_LAMBDA( const int i, const int j ) @@ -331,7 +343,9 @@ class Force> getLinearizedDistanceComponents( x, u, i, j, xi, linear_s, xi_x, xi_y, xi_z ); - const double coeff = c * linear_s * vol( j ); + model.thermalStretch( linear_s, i, j ); + + const double coeff = model.c * linear_s * vol( j ); fx_i = coeff * xi_x / xi; fy_i = coeff * xi_y / xi; fz_i = coeff * xi_z / xi;