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 Oct 31, 2024
1 parent f69e88f commit 2e052bd
Show file tree
Hide file tree
Showing 2 changed files with 103 additions and 43 deletions.
108 changes: 98 additions & 10 deletions src/force/CabanaPD_ForceModels_PMB.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,81 @@ struct ForceModel<PMB, NoFracture, Elastic, TemperatureIndependent>
K = _K;
c = 18.0 * K / ( pi * delta * delta * delta * delta );
}

KOKKOS_INLINE_FUNCTION
bool forceCoeff( const double s, const double vol ) const
{
return c * s * vol;
}

KOKKOS_INLINE_FUNCTION
bool energyCoeff( 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 <>
struct ForceModel<PMB, NoFracture, Plastic, TemperatureIndependent>
: public ForceModel<PMB, NoFracture, Elastic, TemperatureIndependent>
{
using base_type = ForceModel<PMB, NoFracture>;
using base_model = PMB;
using fracture_type = NoFracture;
using plasticity_type = Plastic;
using thermal_type = TemperatureIndependent;

using base_type::c;
using base_type::delta;
using base_type::K;
// FIXME: hardcoded
const double s_Y = 0.0014;

ForceModel(){};
ForceModel( const double delta, const double K )
: base_type( delta, 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 );
}

KOKKOS_INLINE_FUNCTION
bool coeff( const double s, const double vol ) const
{
if ( s < s_Y )
return c * s * vol;
else
return c * s_Y * vol;
}

KOKKOS_INLINE_FUNCTION
bool energyCoeff( 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).
double stretch_term = 0.0;
if ( s < s_Y )
stretch_term = s * s;
else
stretch_term = s_Y * ( 2 * s - s_Y );

return 0.25 * c * stretch_term * xi * vol;
}
};

template <>
Expand Down Expand Up @@ -105,8 +180,7 @@ struct ForceModel<PMB, Fracture, Elastic, TemperatureIndependent>

template <>
struct ForceModel<PMB, Fracture, Plastic, TemperatureIndependent>
: public ForceModel<PMB, Fracture, Elastic, TemperatureIndependent>,
ForceModel<PMB, NoFracture, Plastic, TemperatureIndependent>
: public ForceModel<PMB, Fracture, Elastic, TemperatureIndependent>
{
using base_type = ForceModel<PMB, Fracture>;
using base_model = typename base_type::base_model;
Expand All @@ -120,24 +194,38 @@ struct ForceModel<PMB, Fracture, Plastic, TemperatureIndependent>
using base_type::G0;
using base_type::K;
using base_type::s0;
double alpha = 0.25;

// FIXME: hardcoded
const double s_Y = 0.0014;

ForceModel() {}
ForceModel( const double delta, const double K, const double G0 )
: base_type( delta, K, G0 )
{
}

using base_type::minLocalStretch;
// FIXME: avoiding multiple inheritance..
KOKKOS_INLINE_FUNCTION
bool coeff( const double s, const double vol ) const
{
if ( s < s_Y )
return c * s * vol;
else
return c * s_Y * vol;
}

// Use the dynamic s0 from the minLocalStretch function, reused through the
// neighbor loop.
KOKKOS_INLINE_FUNCTION
bool criticalStretch( const int, const int, const double r, const double xi,
const double s0_i ) const
bool energyCoeff( const double s, const double xi, const double vol ) const
{
double bond_break_coeff = ( 1.0 + s0_i ) * ( 1.0 + s0_i );
return r * r >= bond_break_coeff * xi * xi;
// 0.25 factor is due to 1/2 from outside the integral and 1/2 from
// the integrand (pairwise potential).
double stretch_term = 0.0;
if ( s < s_Y )
stretch_term = s * s;
else
stretch_term = s_Y * ( 2 * s - s_Y );

return 0.25 * c * stretch_term * xi * vol;
}
};

Expand Down
38 changes: 5 additions & 33 deletions src/force/CabanaPD_Force_PMB.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,7 +117,7 @@ class Force<ExecutionSpace, ForceModel<PMB, NoFracture, ModelParams...>>

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;
Expand Down Expand Up @@ -158,7 +158,7 @@ class Force<ExecutionSpace, ForceModel<PMB, NoFracture, ModelParams...>>

// 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.energyCoeff( s, xi, vol( j ) );
W( i ) += w;
Phi += w * vol( i );
};
Expand Down Expand Up @@ -243,19 +243,7 @@ class Force<ExecutionSpace, ForceModel<PMB, Fracture, ModelParams...>>
// Else if statement is only for performance.
else if ( mu( i, n ) > 0 )
{
const double s_Y = 0.0014;
// const double coeff = model.c * s * vol( j );

double coeff = 0;

if ( s < s_Y )
{
coeff = model.c * s * vol( j );
}
else
{
coeff = model.c * s_Y * vol( j );
}
const double coeff = model.forceCoeff( s, vol( j ) );

double muij = mu( i, n );
fx_i = muij * coeff * rx / r;
Expand Down Expand Up @@ -306,23 +294,7 @@ class Force<ExecutionSpace, ForceModel<PMB, Fracture, ModelParams...>>

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
// );

const double s_Y = 0.0014;

if ( s < s_Y )
{
w = mu( i, n ) * 0.25 * model.c * s * s * xi * vol( j );
}
else
{
w = mu( i, n ) * 0.25 * model.c * s_Y * ( 2 * s - s_Y ) *
xi * vol( j );
}

double w = mu( i, n ) * model.energyCoeff( s, xi, vol( j ) );
W( i ) += w;

phi_i += mu( i, n ) * vol( j );
Expand Down Expand Up @@ -391,7 +363,7 @@ class Force<ExecutionSpace, ForceModel<LinearPMB, NoFracture, ModelParams...>>

model.thermalStretch( linear_s, i, j );

const double coeff = model.c * linear_s * vol( j );
const double coeff = model.forceCoeff( linear_s, vol( j ) );
fx_i = coeff * xi_x / xi;
fy_i = coeff * xi_y / xi;
fz_i = coeff * xi_z / xi;
Expand Down

0 comments on commit 2e052bd

Please sign in to comment.