From e5973cd4703727309f094881a677214785e73611 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/2] Remove unused default constructors and set_param --- src/CabanaPD_ForceModels.hpp | 15 ----- src/CabanaPD_Solver.hpp | 1 - src/force/CabanaPD_ContactModels.hpp | 9 --- src/force/CabanaPD_ForceModels_LPS.hpp | 20 +------ src/force/CabanaPD_ForceModels_PMB.hpp | 77 ++------------------------ src/force/CabanaPD_HertzianContact.hpp | 6 -- 6 files changed, 9 insertions(+), 119 deletions(-) diff --git a/src/CabanaPD_ForceModels.hpp b/src/CabanaPD_ForceModels.hpp index 26f600a..78b6a60 100644 --- a/src/CabanaPD_ForceModels.hpp +++ b/src/CabanaPD_ForceModels.hpp @@ -21,12 +21,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 {} @@ -55,12 +52,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. @@ -87,12 +78,6 @@ struct BaseDynamicTemperatureModel BaseDynamicTemperatureModel( const double _delta, const double _kappa, const double _cp, const bool _constant_microconductivity = true ) - { - set_param( _delta, _kappa, _cp, _constant_microconductivity ); - } - - void set_param( const double _delta, const double _kappa, const double _cp, - const bool _constant_microconductivity ) { delta = _delta; kappa = _kappa; diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 6c09158..3aefd4e 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -464,7 +464,6 @@ class SolverNoFracture // Optional modules. std::shared_ptr heat_transfer; std::shared_ptr contact; - contact_model_type contact_model; // Output files. std::string output_file; diff --git a/src/force/CabanaPD_ContactModels.hpp b/src/force/CabanaPD_ContactModels.hpp index 3afe3ae..1f85bef 100644 --- a/src/force/CabanaPD_ContactModels.hpp +++ b/src/force/CabanaPD_ContactModels.hpp @@ -28,7 +28,6 @@ struct ContactModel double delta; double Rc; - ContactModel(){}; // PD horizon // Contact radius ContactModel( const double _delta, const double _Rc ) @@ -51,18 +50,10 @@ struct NormalRepulsionModel : public ContactModel 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 ); diff --git a/src/force/CabanaPD_ForceModels_LPS.hpp b/src/force/CabanaPD_ForceModels_LPS.hpp index 2317b26..45f2d1f 100644 --- a/src/force/CabanaPD_ForceModels_LPS.hpp +++ b/src/force/CabanaPD_ForceModels_LPS.hpp @@ -40,16 +40,9 @@ struct ForceModel : public BaseForceModel 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; } @@ -85,15 +78,8 @@ struct 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 = Kokkos::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 3b6180b..9e278c1 100644 --- a/src/force/CabanaPD_ForceModels_PMB.hpp +++ b/src/force/CabanaPD_ForceModels_PMB.hpp @@ -34,23 +34,10 @@ struct ForceModel double c; double K; - 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 ); } }; @@ -71,24 +58,10 @@ struct ForceModel double s0; double bond_break_coeff; - 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 = Kokkos::sqrt( 5.0 * G0 / 9.0 / K / delta ); bond_break_coeff = ( 1.0 + s0 ) * ( 1.0 + s0 ); } @@ -167,14 +140,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 @@ -301,17 +258,6 @@ struct ForceModel , base_temperature_type( _delta, _kappa, _cp, _constant_microconductivity ) { - set_param( _delta, _K, _kappa, _cp, _alpha, _temp0, - _constant_microconductivity ); - } - - void set_param( const double _delta, const double _K, const double _kappa, - const double _cp, const double _alpha, const double _temp0, - const bool _constant_microconductivity ) - { - base_type::set_param( _delta, _K, _alpha, _temp0 ); - base_temperature_type::set_param( _delta, _kappa, _cp, - _constant_microconductivity ); } }; @@ -364,20 +310,9 @@ struct ForceModel const double _temp0 = 0.0, const bool _constant_microconductivity = true ) : base_type( _delta, _K, _G0, _temp, _alpha, _temp0 ) + , base_temperature_type( _delta, _kappa, _cp, + _constant_microconductivity ) { - set_param( _delta, _K, _G0, _kappa, _cp, _alpha, _temp0 ); - base_temperature_type::set_param( _delta, _kappa, _cp, - _constant_microconductivity ); - } - - void set_param( const double _delta, const double _K, const double _G0, - const double _kappa, const double _cp, const double _alpha, - const double _temp0, - const bool _constant_microconductivity ) - { - base_type::set_param( _delta, _K, _G0, _alpha, _temp0 ); - base_temperature_type::set_param( _delta, _kappa, _cp, - _constant_microconductivity ); } }; diff --git a/src/force/CabanaPD_HertzianContact.hpp b/src/force/CabanaPD_HertzianContact.hpp index b99aac9..0ad012d 100644 --- a/src/force/CabanaPD_HertzianContact.hpp +++ b/src/force/CabanaPD_HertzianContact.hpp @@ -39,12 +39,6 @@ struct HertzianModel : public ContactModel HertzianModel( const double _Rc, const double _radius, const double _nu, const double _E, const double _e ) : ContactModel( 1.0, _Rc ) - { - set_param( _radius, _nu, _E, _e ); - } - - void set_param( const double _radius, const double _nu, const double _E, - const double _e ) { nu = _nu; radius = _radius; From 22ce38389ec283e782476602eaea435f4214bee3 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Thu, 9 Jan 2025 12:35:39 -0500 Subject: [PATCH 2/2] Add model interfaces for force/energy calculations --- src/force/CabanaPD_ContactModels.hpp | 9 +++ src/force/CabanaPD_ForceModels_LPS.hpp | 38 +++++++++- src/force/CabanaPD_ForceModels_PMB.hpp | 14 ++++ src/force/CabanaPD_Force_Contact.hpp | 13 +--- src/force/CabanaPD_Force_HertzianContact.hpp | 59 ++++----------- src/force/CabanaPD_Force_LPS.hpp | 79 ++++++-------------- src/force/CabanaPD_Force_PMB.hpp | 19 ++--- src/force/CabanaPD_HertzianContact.hpp | 35 +++++++++ unit_test/tstForce.hpp | 15 ++-- 9 files changed, 149 insertions(+), 132 deletions(-) diff --git a/src/force/CabanaPD_ContactModels.hpp b/src/force/CabanaPD_ContactModels.hpp index 1f85bef..aff456f 100644 --- a/src/force/CabanaPD_ContactModels.hpp +++ b/src/force/CabanaPD_ContactModels.hpp @@ -58,6 +58,15 @@ struct NormalRepulsionModel : public ContactModel // This could inherit from PMB (same c) c = 18.0 * K / ( pi * delta * delta * delta * delta ); } + + KOKKOS_INLINE_FUNCTION + auto forceCoeff( const double r, const double vol ) const + { + // Contact "stretch" + const double sc = ( r - Rc ) / delta; + // Normal repulsion uses a 15 factor compared to the PMB force + return 15.0 * c * sc * vol; + } }; template <> diff --git a/src/force/CabanaPD_ForceModels_LPS.hpp b/src/force/CabanaPD_ForceModels_LPS.hpp index 45f2d1f..f8538de 100644 --- a/src/force/CabanaPD_ForceModels_LPS.hpp +++ b/src/force/CabanaPD_ForceModels_LPS.hpp @@ -47,13 +47,49 @@ struct ForceModel : public BaseForceModel s_coeff = 15.0 * G; } - KOKKOS_INLINE_FUNCTION double influence_function( double xi ) const + KOKKOS_INLINE_FUNCTION double influenceFunction( double xi ) const { if ( influence_type == 1 ) return 1.0 / xi; else return 1.0; } + + KOKKOS_INLINE_FUNCTION auto weightedVolume( const double xi, + const double vol ) const + { + return influenceFunction( xi ) * xi * xi * vol; + } + + KOKKOS_INLINE_FUNCTION auto dilatation( const double s, const double xi, + const double vol, + const double m_i ) const + { + double theta_i = influenceFunction( xi ) * s * xi * xi * vol; + return 3.0 * theta_i / m_i; + } + + KOKKOS_INLINE_FUNCTION auto forceCoeff( const double s, const double xi, + const double vol, const double m_i, + const double m_j, + const double theta_i, + const double theta_j ) const + { + return ( theta_coeff * ( theta_i / m_i + theta_j / m_j ) + + s_coeff * s * ( 1.0 / m_i + 1.0 / m_j ) ) * + influenceFunction( xi ) * xi * vol; + } + + KOKKOS_INLINE_FUNCTION + auto energy( const double s, const double xi, const double vol, + const double m_i, const double theta_i, + const double num_bonds ) const + { + return 1.0 / num_bonds * 0.5 * theta_coeff / 3.0 * + ( theta_i * theta_i ) + + 0.5 * ( s_coeff / m_i ) * influenceFunction( xi ) * s * s * xi * + xi * vol; + } }; template <> diff --git a/src/force/CabanaPD_ForceModels_PMB.hpp b/src/force/CabanaPD_ForceModels_PMB.hpp index 9e278c1..aba3860 100644 --- a/src/force/CabanaPD_ForceModels_PMB.hpp +++ b/src/force/CabanaPD_ForceModels_PMB.hpp @@ -40,6 +40,20 @@ struct ForceModel { c = 18.0 * K / ( pi * delta * delta * delta * delta ); } + + KOKKOS_INLINE_FUNCTION + auto forceCoeff( const double s, const double vol ) const + { + return c * s * vol; + } + + KOKKOS_INLINE_FUNCTION + auto energy( const double s, const double xi, const double vol ) const + { + // 0.25 factor is due to 1/2 from outside the integral and 1/2 from + // the integrand (pairwise potential). + return 0.25 * c * s * s * xi * vol; + } }; template <> diff --git a/src/force/CabanaPD_Force_Contact.hpp b/src/force/CabanaPD_Force_Contact.hpp index ba14363..9f299f5 100644 --- a/src/force/CabanaPD_Force_Contact.hpp +++ b/src/force/CabanaPD_Force_Contact.hpp @@ -69,16 +69,15 @@ class Force const ParticleType& particles, ParallelType& neigh_op_tag ) { - auto delta = _model.delta; - auto Rc = _model.Rc; - auto c = _model.c; + auto model = _model; const auto vol = particles.sliceVolume(); const auto y = particles.sliceCurrentPosition(); const int n_frozen = particles.frozenOffset(); const int n_local = particles.localOffset(); _neigh_timer.start(); - _neigh_list.build( y, n_frozen, n_local, Rc, 1.0, mesh_min, mesh_max ); + _neigh_list.build( y, n_frozen, n_local, model.Rc, 1.0, mesh_min, + mesh_max ); _neigh_timer.stop(); auto contact_full = KOKKOS_LAMBDA( const int i, const int j ) @@ -91,11 +90,7 @@ class Force 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 ); + const double coeff = model.forceCoeff( r, vol( j ) ); fcx_i = coeff * rx / r; fcy_i = coeff * ry / r; fcz_i = coeff * rz / r; diff --git a/src/force/CabanaPD_Force_HertzianContact.hpp b/src/force/CabanaPD_Force_HertzianContact.hpp index 688636b..1cdde1d 100644 --- a/src/force/CabanaPD_Force_HertzianContact.hpp +++ b/src/force/CabanaPD_Force_HertzianContact.hpp @@ -54,75 +54,44 @@ class Force const ParticleType& particles, ParallelType& neigh_op_tag ) { - auto Rc = _model.Rc; - auto radius = _model.radius; - auto Es = _model.Es; - auto Rs = _model.Rs; - auto beta = _model.beta; const int n_frozen = particles.frozenOffset(); const int n_local = particles.localOffset(); - const double coeff_h_n = 4.0 / 3.0 * Es * Kokkos::sqrt( Rs ); - const double coeff_h_d = -2.0 * Kokkos::sqrt( 5.0 / 6.0 ) * beta; - + auto model = _model; const auto vol = particles.sliceVolume(); const auto rho = particles.sliceDensity(); const auto y = particles.sliceCurrentPosition(); const auto vel = particles.sliceVelocity(); _neigh_timer.start(); - _neigh_list.build( y, n_frozen, n_local, Rc, 1.0, mesh_min, mesh_max ); + _neigh_list.build( y, n_frozen, n_local, model.Rc, 1.0, mesh_min, + mesh_max ); _neigh_timer.stop(); 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 "overlap" - const double delta_n = ( r - 2.0 * radius ); - - // Hertz normal force coefficient - double coeff = 0.0; - double Sn = 0.0; - if ( delta_n < 0.0 ) - { - coeff = Kokkos::min( - 0.0, -coeff_h_n * - Kokkos::pow( Kokkos::abs( delta_n ), 3.0 / 2.0 ) ); - Sn = 2.0 * Es * Kokkos::sqrt( Rs * Kokkos::abs( delta_n ) ); - } + const double delta_n = ( r - 2.0 * model.radius ); - coeff /= vol( i ); - - 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; + const double coeff_normal = + model.normalForceCoeff( delta_n, vol( i ) ); + fc( i, 0 ) += coeff_normal * rx / r; + fc( i, 1 ) += coeff_normal * ry / r; + fc( i, 2 ) += coeff_normal * rz / r; // Hertz normal force damping component double vx, vy, vz, vn; getRelativeNormalVelocityComponents( vel, i, j, rx, ry, rz, r, vx, vy, vz, vn ); - double ms = ( rho( i ) * vol( i ) ) / 2.0; - double fnd = coeff_h_d * Kokkos::sqrt( Sn * ms ) * vn / vol( i ); - - fcx_i = fnd * rx / r; - fcy_i = fnd * ry / r; - fcz_i = fnd * rz / r; - - fc( i, 0 ) += fcx_i; - fc( i, 1 ) += fcy_i; - fc( i, 2 ) += fcz_i; + const double coeff_damp = + model.dampingForceCoeff( delta_n, vn, vol( i ), rho( i ) ); + fc( i, 0 ) += coeff_damp * rx / r; + fc( i, 1 ) += coeff_damp * ry / r; + fc( i, 2 ) += coeff_damp * rz / r; }; _timer.start(); diff --git a/src/force/CabanaPD_Force_LPS.hpp b/src/force/CabanaPD_Force_LPS.hpp index 8262e62..a3f5747 100644 --- a/src/force/CabanaPD_Force_LPS.hpp +++ b/src/force/CabanaPD_Force_LPS.hpp @@ -114,8 +114,7 @@ class Force> // Get the reference positions and displacements. double xi, r, s; getDistance( x, u, i, j, xi, r, s ); - double m_j = model.influence_function( xi ) * xi * xi * vol( j ); - m( i ) += m_j; + m( i ) += model.weightedVolume( xi, vol( j ) ); }; Kokkos::RangePolicy policy( particles.frozenOffset(), @@ -146,9 +145,7 @@ class Force> // Get the bond distance, displacement, and stretch. double xi, r, s; getDistance( x, u, i, j, xi, r, s ); - double theta_i = - model.influence_function( xi ) * s * xi * xi * vol( j ); - theta( i ) += 3.0 * theta_i / m( i ); + theta( i ) += model.dilatation( s, xi, vol( j ), m( i ) ); }; Kokkos::RangePolicy policy( particles.frozenOffset(), @@ -168,10 +165,7 @@ class Force> { _timer.start(); - auto theta_coeff = _model.theta_coeff; - auto s_coeff = _model.s_coeff; auto model = _model; - const auto vol = particles.sliceVolume(); auto theta = particles.sliceDilatation(); auto m = particles.sliceWeightedVolume(); @@ -184,10 +178,9 @@ class Force> double xi, r, s; double rx, ry, rz; getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz ); - const double coeff = - ( theta_coeff * ( theta( i ) / m( i ) + theta( j ) / m( j ) ) + - s_coeff * s * ( 1.0 / m( i ) + 1.0 / m( j ) ) ) * - model.influence_function( xi ) * xi * vol( j ); + + const double coeff = model.forceCoeff( + s, xi, vol( j ), m( i ), m( j ), theta( i ), theta( j ) ); fx_i = coeff * rx / r; fy_i = coeff * ry / r; fz_i = coeff * rz / r; @@ -214,8 +207,6 @@ class Force> { _energy_timer.start(); - auto theta_coeff = _model.theta_coeff; - auto s_coeff = _model.s_coeff; auto model = _model; auto neigh_list = _neigh_list; @@ -233,11 +224,8 @@ class Force> Cabana::NeighborList::numNeighbor( neigh_list, i ) ); - double w = ( 1.0 / num_neighbors ) * 0.5 * theta_coeff / 3.0 * - ( theta( i ) * theta( i ) ) + - 0.5 * ( s_coeff / m( i ) ) * - model.influence_function( xi ) * s * s * xi * xi * - vol( j ); + double w = model.energy( s, xi, vol( j ), m( i ), theta( i ), + num_neighbors ); W( i ) += w; Phi += w * vol( i ); }; @@ -314,9 +302,7 @@ class Force> double xi, r, s; getDistance( x, u, i, j, xi, r, s ); // mu is included to account for bond breaking. - double m_j = mu( i, n ) * model.influence_function( xi ) * xi * - xi * vol( j ); - m( i ) += m_j; + m( i ) += mu( i, n ) * model.weightedVolume( xi, vol( j ) ); } }; @@ -356,14 +342,14 @@ class Force> // Get the bond distance, displacement, and stretch. double xi, r, s; getDistance( x, u, i, j, xi, r, s ); - // mu is included to account for bond breaking. - double theta_i = mu( i, n ) * model.influence_function( xi ) * - s * xi * xi * vol( j ); + // Check if all bonds are broken (m=0) to avoid dividing by // zero. Alternatively, one could check if this bond mu(i,n) is // broken, because m=0 only occurs when all bonds are broken. + // mu is still included to account for individual bond breaking. if ( m( i ) > 0 ) - theta( i ) += 3.0 * theta_i / m( i ); + theta( i ) += mu( i, n ) * + model.dilatation( s, xi, vol( j ), m( i ) ); } }; @@ -384,8 +370,6 @@ class Force> _timer.start(); auto break_coeff = _model.bond_break_coeff; - auto theta_coeff = _model.theta_coeff; - auto s_coeff = _model.s_coeff; auto model = _model; auto neigh_list = _neigh_list; @@ -425,10 +409,8 @@ class Force> else if ( mu( i, n ) > 0 ) { const double coeff = - ( theta_coeff * - ( theta( i ) / m( i ) + theta( j ) / m( j ) ) + - s_coeff * s * ( 1.0 / m( i ) + 1.0 / m( j ) ) ) * - model.influence_function( xi ) * xi * vol( j ); + model.forceCoeff( s, xi, vol( j ), m( i ), m( j ), + theta( i ), theta( j ) ); double muij = mu( i, n ); fx_i = muij * coeff * rx / r; fy_i = muij * coeff * ry / r; @@ -457,8 +439,6 @@ class Force> { _energy_timer.start(); - auto theta_coeff = _model.theta_coeff; - auto s_coeff = _model.s_coeff; auto model = _model; auto neigh_list = _neigh_list; @@ -486,12 +466,8 @@ class Force> double xi, r, s; getDistance( x, u, i, j, xi, r, s ); - double w = - mu( i, n ) * ( ( 1.0 / num_bonds ) * 0.5 * theta_coeff / - 3.0 * ( theta( i ) * theta( i ) ) + - 0.5 * ( s_coeff / m( i ) ) * - model.influence_function( xi ) * s * s * - xi * xi * vol( j ) ); + double w = mu( i, n ) * model.energy( s, xi, vol( j ), m( i ), + theta( i ), num_bonds ); W( i ) += w; phi_i += mu( i, n ) * vol( j ); @@ -546,12 +522,9 @@ class Force> { _timer.start(); - auto theta_coeff = _model.theta_coeff; - auto s_coeff = _model.s_coeff; auto model = _model; - const auto vol = particles.sliceVolume(); - auto linear_theta = particles.sliceDilatation(); + auto theta = particles.sliceDilatation(); // Using weighted volume from base LPS class. auto m = particles.sliceWeightedVolume(); @@ -568,10 +541,8 @@ class Force> xi_y, xi_z ); const double coeff = - ( theta_coeff * ( linear_theta( i ) / m( i ) + - linear_theta( j ) / m( j ) ) + - s_coeff * linear_s * ( 1.0 / m( i ) + 1.0 / m( j ) ) ) * - model.influence_function( xi ) * xi * vol( j ); + model.forceCoeff( linear_s, xi, vol( j ), m( i ), m( j ), + theta( i ), theta( j ) ); fx_i = coeff * xi_x / xi; fy_i = coeff * xi_y / xi; fz_i = coeff * xi_z / xi; @@ -598,20 +569,17 @@ class Force> { _energy_timer.start(); - auto theta_coeff = _model.theta_coeff; - auto s_coeff = _model.s_coeff; auto model = _model; auto neigh_list = _neigh_list; const auto vol = particles.sliceVolume(); - const auto linear_theta = particles.sliceDilatation(); + const auto theta = particles.sliceDilatation(); auto m = particles.sliceWeightedVolume(); auto energy_full = KOKKOS_LAMBDA( const int i, const int j, double& Phi ) { // Do we need to recompute linear_theta_i? - double xi, linear_s; getLinearizedDistance( x, u, i, j, xi, linear_s ); @@ -619,11 +587,8 @@ class Force> Cabana::NeighborList::numNeighbor( neigh_list, i ) ); - double w = ( 1.0 / num_neighbors ) * 0.5 * theta_coeff / 3.0 * - ( linear_theta( i ) * linear_theta( i ) ) + - 0.5 * ( s_coeff / m( i ) ) * - model.influence_function( xi ) * linear_s * - linear_s * xi * xi * vol( j ); + double w = model.energy( linear_s, xi, vol( j ), m( i ), theta( i ), + num_neighbors ); W( i ) += w; Phi += w * vol( i ); }; diff --git a/src/force/CabanaPD_Force_PMB.hpp b/src/force/CabanaPD_Force_PMB.hpp index 8a09963..4079851 100644 --- a/src/force/CabanaPD_Force_PMB.hpp +++ b/src/force/CabanaPD_Force_PMB.hpp @@ -121,7 +121,7 @@ class Force> model.thermalStretch( s, i, j ); - const double coeff = model.c * s * vol( j ); + const double coeff = model.forceCoeff( s, vol( j ) ); fx_i = coeff * rx / r; fy_i = coeff * ry / r; fz_i = coeff * rz / r; @@ -160,9 +160,7 @@ class Force> model.thermalStretch( s, i, j ); - // 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 = model.energy( s, xi, vol( j ) ); W( i ) += w; Phi += w * vol( i ); }; @@ -253,7 +251,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.forceCoeff( s, vol( j ) ); + double muij = mu( i, n ); fx_i = muij * coeff * rx / r; fy_i = muij * coeff * ry / r; @@ -304,9 +303,7 @@ class Force> model.thermalStretch( s, i, j ); - // 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 ) * model.energy( s, xi, vol( j ) ); W( i ) += w; phi_i += mu( i, n ) * vol( j ); @@ -381,7 +378,7 @@ class Force 0.0 && xi < model.delta + 1e-14 ) weighted_volume += - model.influence_function( xi ) * xi * xi * vol; + model.influenceFunction( xi ) * xi * xi * vol; } return weighted_volume; } @@ -198,8 +198,7 @@ double computeReferenceDilatation( ModelType model, const int m, if ( xi > 0.0 && xi < model.delta + 1e-14 ) theta += 3.0 / weighted_volume * - model.influence_function( xi ) * s0 * xi * xi * - vol; + model.influenceFunction( xi ) * s0 * xi * xi * vol; } return theta; } @@ -231,7 +230,7 @@ double computeReferenceDilatation( ModelType model, const int m, if ( xi > 0.0 && xi < model.delta + 1e-14 ) theta += 3.0 / weighted_volume * - model.influence_function( xi ) * s * xi * xi * vol; + model.influenceFunction( xi ) * s * xi * xi * vol; } return theta; } @@ -287,8 +286,8 @@ double computeReferenceStrainEnergyDensity( W += ( 1.0 / num_neighbors ) * 0.5 * model.theta_coeff / 3.0 * ( theta * theta ) + 0.5 * ( model.s_coeff / weighted_volume ) * - model.influence_function( xi ) * s0 * s0 * xi * - xi * vol; + model.influenceFunction( xi ) * s0 * s0 * xi * xi * + vol; } } return W; @@ -335,7 +334,7 @@ double computeReferenceStrainEnergyDensity( W += ( 1.0 / num_neighbors ) * 0.5 * model.theta_coeff / 3.0 * ( theta_i * theta_j ) + 0.5 * ( model.s_coeff / weighted_volume ) * - model.influence_function( xi ) * s * s * xi * xi * + model.influenceFunction( xi ) * s * s * xi * xi * vol; } } @@ -384,7 +383,7 @@ double computeReferenceForceX( model.s_coeff * s * ( 1.0 / weighted_volume + 1.0 / weighted_volume ) ) * - model.influence_function( xi ) * xi * vol * rx / r; + model.influenceFunction( xi ) * xi * vol * rx / r; } } return fx;