Skip to content

Commit

Permalink
Add model interfaces for force/energy calculations
Browse files Browse the repository at this point in the history
  • Loading branch information
@streeve
streeve committed Jan 9, 2025
1 parent e5973cd commit 22ce383
Show file tree
Hide file tree
Showing 9 changed files with 149 additions and 132 deletions.
9 changes: 9 additions & 0 deletions src/force/CabanaPD_ContactModels.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 <>
Expand Down
38 changes: 37 additions & 1 deletion src/force/CabanaPD_ForceModels_LPS.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,13 +47,49 @@ struct ForceModel<LPS, Elastic, NoFracture> : 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 <>
Expand Down
14 changes: 14 additions & 0 deletions src/force/CabanaPD_ForceModels_PMB.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,20 @@ struct ForceModel<PMB, Elastic, NoFracture, TemperatureIndependent>
{
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 <>
Expand Down
13 changes: 4 additions & 9 deletions src/force/CabanaPD_Force_Contact.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,16 +69,15 @@ class Force<MemorySpace, NormalRepulsionModel>
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 )
Expand All @@ -91,11 +90,7 @@ class Force<MemorySpace, NormalRepulsionModel>
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;
Expand Down
59 changes: 14 additions & 45 deletions src/force/CabanaPD_Force_HertzianContact.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,75 +54,44 @@ class Force<MemorySpace, HertzianModel>
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();
Expand Down
Loading

0 comments on commit 22ce383

Please sign in to comment.