tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Minimizing changes to pair granular

Browse Source
This commit is contained in:
jtclemm
2024-12-18 17:44:12 -07:00
parent d764c367c7
commit 844d3a6cab
6 changed files with 231 additions and 295 deletions
Download Patch File Download Diff File Expand all files Collapse all files

306
src/GRANULAR/fix_granular_mdr.cpp
View File

@ -54,7 +54,15 @@ enum {MEAN_SURF_DISP, RADIUS_UPDATE};
FixGranularMDR::FixGranularMDR(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
comm_forward = 20; // value needs to match number of values you communicate
comm_forward = 20; // value needs to match number of values you communicate
id_fix = nullptr;
}
/* ---------------------------------------------------------------------- */
FixGranularMDR::~FixGranularMDR()
{
if (id_fix) modify->delete_fix(id_fix);
}
/* ---------------------------------------------------------------------- */
@ -68,13 +76,16 @@ int FixGranularMDR::setmask()
/* ---------------------------------------------------------------------- */
void FixGranularMDR::setup_pre_force(int /*vflag*/)
void FixGranularMDR::post_constructor()
{
int tmp1, tmp2;
id_fix = utils::strdup("MDR_PARTICLE_HISTORY_VARIABLES");
modify->add_fix(fmt::format("{} all property/atom d_Ro d_Vcaps d_Vgeo d_Velas d_eps_bar d_dRnumerator d_dRdenominator d_Acon0 d_Acon1 d_Atot d_Atot_sum d_ddelta_bar d_psi d_psi_b d_history_setup_flag d_sigmaxx d_sigmayy d_sigmazz d_contacts d_adhesive_length ghost yes", id_fix));
int index_Ro = atom->find_custom("Ro",tmp1,tmp2);
int index_Vcaps = atom->find_custom("Vcaps",tmp1,tmp2);
int index_Vgeo = atom->find_custom("Vgeo",tmp1,tmp2);
int index_Velas = atom->find_custom("Velas",tmp1,tmp2);
int index_Vcaps = atom->find_custom("Vcaps",tmp1,tmp2);
int index_eps_bar = atom->find_custom("eps_bar",tmp1,tmp2);
int index_dRnumerator = atom->find_custom("dRnumerator",tmp1,tmp2);
int index_dRdenominator = atom->find_custom("dRdenominator",tmp1,tmp2);
@ -85,12 +96,13 @@ void FixGranularMDR::setup_pre_force(int /*vflag*/)
int index_ddelta_bar = atom->find_custom("ddelta_bar",tmp1,tmp2);
int index_psi = atom->find_custom("psi",tmp1,tmp2);
int index_psi_b = atom->find_custom("psi_b",tmp1,tmp2);
int index_history_setup_flag = atom->find_custom("history_setup_flag",tmp1,tmp2);
int index_sigmaxx = atom->find_custom("sigmaxx",tmp1,tmp2);
int index_sigmayy = atom->find_custom("sigmayy",tmp1,tmp2);
int index_sigmazz = atom->find_custom("sigmazz",tmp1,tmp2);
int index_history_setup_flag = atom->find_custom("history_setup_flag",tmp1,tmp2);
int index_contacts = atom->find_custom("contacts",tmp1,tmp2);
int index_adhesive_length = atom->find_custom("adhesive_length",tmp1,tmp2);
Ro = atom->dvector[index_Ro];
Vgeo = atom->dvector[index_Vgeo];
Velas = atom->dvector[index_Velas];
@ -111,7 +123,12 @@ void FixGranularMDR::setup_pre_force(int /*vflag*/)
history_setup_flag = atom->dvector[index_history_setup_flag];
contacts = atom->dvector[index_contacts];
adhesive_length = atom->dvector[index_adhesive_length];
}
/* ---------------------------------------------------------------------- */
void FixGranularMDR::setup_pre_force(int /*vflag*/)
{
pre_force(0);
}
@ -119,48 +136,6 @@ void FixGranularMDR::setup_pre_force(int /*vflag*/)
void FixGranularMDR::setup(int /*vflag*/)
{
int tmp1, tmp2;
int index_Ro = atom->find_custom("Ro",tmp1,tmp2);
int index_Vgeo = atom->find_custom("Vgeo",tmp1,tmp2);
int index_Velas = atom->find_custom("Velas",tmp1,tmp2);
int index_Vcaps = atom->find_custom("Vcaps",tmp1,tmp2);
int index_eps_bar = atom->find_custom("eps_bar",tmp1,tmp2);
int index_dRnumerator = atom->find_custom("dRnumerator",tmp1,tmp2);
int index_dRdenominator = atom->find_custom("dRdenominator",tmp1,tmp2);
int index_Acon0 = atom->find_custom("Acon0",tmp1,tmp2);
int index_Acon1 = atom->find_custom("Acon1",tmp1,tmp2);
int index_Atot = atom->find_custom("Atot",tmp1,tmp2);
int index_Atot_sum = atom->find_custom("Atot_sum",tmp1,tmp2);
int index_ddelta_bar = atom->find_custom("ddelta_bar",tmp1,tmp2);
int index_psi = atom->find_custom("psi",tmp1,tmp2);
int index_psi_b = atom->find_custom("psi_b",tmp1,tmp2);
int index_sigmaxx = atom->find_custom("sigmaxx",tmp1,tmp2);
int index_sigmayy = atom->find_custom("sigmayy",tmp1,tmp2);
int index_sigmazz = atom->find_custom("sigmazz",tmp1,tmp2);
int index_history_setup_flag = atom->find_custom("history_setup_flag",tmp1,tmp2);
int index_contacts = atom->find_custom("contacts",tmp1,tmp2);
int index_adhesive_length = atom->find_custom("adhesive_length",tmp1,tmp2);
Ro = atom->dvector[index_Ro];
Vgeo = atom->dvector[index_Vgeo];
Velas = atom->dvector[index_Velas];
Vcaps = atom->dvector[index_Vcaps];
eps_bar = atom->dvector[index_eps_bar];
dRnumerator = atom->dvector[index_dRnumerator];
dRdenominator = atom->dvector[index_dRdenominator];
Acon0 = atom->dvector[index_Acon0];
Acon1 = atom->dvector[index_Acon1];
Atot = atom->dvector[index_Atot];
Atot_sum = atom->dvector[index_Atot_sum];
ddelta_bar = atom->dvector[index_ddelta_bar];
psi = atom->dvector[index_psi];
psi_b = atom->dvector[index_psi_b];
sigmaxx = atom->dvector[index_sigmaxx];
sigmayy = atom->dvector[index_sigmayy];
sigmazz = atom->dvector[index_sigmazz];
history_setup_flag = atom->dvector[index_history_setup_flag];
contacts = atom->dvector[index_contacts];
adhesive_length = atom->dvector[index_adhesive_length];
end_of_step();
}
@ -381,121 +356,121 @@ void FixGranularMDR::calculate_contact_penalty()
double radi = radius[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (int jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
double radj = radius[j];
const double delx_ij = x[j][0] - xtmp;
const double dely_ij = x[j][1] - ytmp;
const double delz_ij = x[j][2] - ztmp;
const double rsq_ij = delx_ij*delx_ij + dely_ij*dely_ij + delz_ij*delz_ij;
const double r_ij = sqrt(rsq_ij);
const double rinv_ij = 1.0/r_ij;
const double radsum_ij = radi + radj;
const double deltan_ij = radsum_ij - r_ij;
if (deltan_ij >= 0.0) {
for (int kk = jj; kk < jnum; kk++) {
k = jlist[kk];
k &= NEIGHMASK;
ktype = type[k];
if (kk != jj) {
const double delx_ik = x[k][0] - xtmp;
const double dely_ik = x[k][1] - ytmp;
const double delz_ik = x[k][2] - ztmp;
const double rsq_ik = delx_ik*delx_ik + dely_ik*dely_ik + delz_ik*delz_ik;
const double r_ik = sqrt(rsq_ik);
const double rinv_ik = 1.0/r_ik;
const double radk = radius[k];
const double radsum_ik = radi + radk;
const double deltan_ik = radsum_ik - r_ik;
const double delx_jk = x[k][0] - x[j][0];
const double dely_jk = x[k][1] - x[j][1];
const double delz_jk = x[k][2] - x[j][2];
const double rsq_jk = delx_jk*delx_jk + dely_jk*dely_jk + delz_jk*delz_jk;
const double r_jk = sqrt(rsq_jk);
const double rinv_jk = 1.0/r_jk;
const double radsum_jk = radj + radk;
const double deltan_jk = radsum_jk - r_jk;
if (deltan_ik >= 0.0 && deltan_jk >= 0.0) {
for (int jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
double radj = radius[j];
const double delx_ij = x[j][0] - xtmp;
const double dely_ij = x[j][1] - ytmp;
const double delz_ij = x[j][2] - ztmp;
const double rsq_ij = delx_ij*delx_ij + dely_ij*dely_ij + delz_ij*delz_ij;
const double r_ij = sqrt(rsq_ij);
const double rinv_ij = 1.0/r_ij;
const double radsum_ij = radi + radj;
const double deltan_ij = radsum_ij - r_ij;
if (deltan_ij >= 0.0) {
for (int kk = jj; kk < jnum; kk++) {
k = jlist[kk];
k &= NEIGHMASK;
ktype = type[k];
if (kk != jj) {
const double delx_ik = x[k][0] - xtmp;
const double dely_ik = x[k][1] - ytmp;
const double delz_ik = x[k][2] - ztmp;
const double rsq_ik = delx_ik*delx_ik + dely_ik*dely_ik + delz_ik*delz_ik;
const double r_ik = sqrt(rsq_ik);
const double rinv_ik = 1.0/r_ik;
const double radk = radius[k];
const double radsum_ik = radi + radk;
const double deltan_ik = radsum_ik - r_ik;
const double delx_jk = x[k][0] - x[j][0];
const double dely_jk = x[k][1] - x[j][1];
const double delz_jk = x[k][2] - x[j][2];
const double rsq_jk = delx_jk*delx_jk + dely_jk*dely_jk + delz_jk*delz_jk;
const double r_jk = sqrt(rsq_jk);
const double rinv_jk = 1.0/r_jk;
const double radsum_jk = radj + radk;
const double deltan_jk = radsum_jk - r_jk;
if (deltan_ik >= 0.0 && deltan_jk >= 0.0) {
// pull ij history
history_ij = &allhistory[size_history * jj];
double * pij = &history_ij[22]; // penalty for contact i and j
// pull ij history
history_ij = &allhistory[size_history * jj];
double * pij = &history_ij[22]; // penalty for contact i and j
// pull ik history
history_ik = &allhistory[size_history * kk];
double * pik = &history_ik[22]; // penalty for contact i and k
// pull ik history
history_ik = &allhistory[size_history * kk];
double * pik = &history_ik[22]; // penalty for contact i and k
// we don't know if who owns the contact ahead of time, k might be in j's neigbor list or vice versa,
// so we need to manually search to figure out the owner check if k is in the neighbor list of j
double * pjk = NULL;
int * const jklist = firstneigh[j];
const int jknum = numneigh[j];
for (int jk = 0; jk < jknum; jk++) {
const int kneigh = jklist[jk] & NEIGHMASK;
if (k == kneigh) {
allhistory_j = firsthistory[j];
history_jk = &allhistory_j[size_history * jk];
pjk = &history_jk[22]; // penalty for contact j and k
// we don't know if who owns the contact ahead of time, k might be in j's neigbor list or vice versa,
// so we need to manually search to figure out the owner check if k is in the neighbor list of j
double * pjk = NULL;
int * const jklist = firstneigh[j];
const int jknum = numneigh[j];
for (int jk = 0; jk < jknum; jk++) {
const int kneigh = jklist[jk] & NEIGHMASK;
if (k == kneigh) {
allhistory_j = firsthistory[j];
history_jk = &allhistory_j[size_history * jk];
pjk = &history_jk[22]; // penalty for contact j and k
break;
}
}
// check if j is in the neighbor list of k
if (pjk == NULL) {
int * const kjlist = firstneigh[k];
const int kjnum = numneigh[k];
for (int kj = 0; kj < kjnum; kj++) {
const int jneigh = kjlist[kj] & NEIGHMASK;
if (j == jneigh) {
allhistory_k = firsthistory[k];
history_kj = &allhistory_k[size_history * kj];
pjk = &history_kj[22]; // penalty for contact j and k
break;
}
}
}
// check if j is in the neighbor list of k
if (pjk == NULL) {
int * const kjlist = firstneigh[k];
const int kjnum = numneigh[k];
for (int kj = 0; kj < kjnum; kj++) {
const int jneigh = kjlist[kj] & NEIGHMASK;
if (j == jneigh) {
allhistory_k = firsthistory[k];
history_kj = &allhistory_k[size_history * kj];
pjk = &history_kj[22]; // penalty for contact j and k
break;
}
}
}
std::vector<double> distances = {r_ij,r_ik,r_jk};
auto maxElement = std::max_element(distances.begin(), distances.end());
double maxValue = *maxElement;
int maxIndex = std::distance(distances.begin(), maxElement);
if (maxIndex == 0) { // the central particle is k
const double enx_ki = -delx_ik * rinv_ik;
const double eny_ki = -dely_ik * rinv_ik;
const double enz_ki = -delz_ik * rinv_ik;
const double enx_kj = -delx_jk * rinv_jk;
const double eny_kj = -dely_jk * rinv_jk;
const double enz_kj = -delz_jk * rinv_jk;
const double alpha = std::acos(enx_ki*enx_kj + eny_ki*eny_kj + enz_ki*enz_kj);
pij[0] += 1.0/( 1.0 + std::exp(-50.0*(alpha/MY_PI - 1.0/2.0)) );
} else if (maxIndex == 1) { // the central particle is j
const double enx_ji = -delx_ij * rinv_ij;
const double eny_ji = -dely_ij * rinv_ij;
const double enz_ji = -delz_ij * rinv_ij;
const double enx_jk = delx_jk * rinv_jk;
const double eny_jk = dely_jk * rinv_jk;
const double enz_jk = delz_jk * rinv_jk;
const double alpha = std::acos(enx_ji*enx_jk + eny_ji*eny_jk + enz_ji*enz_jk);
pik[0] += 1.0/( 1.0 + std::exp(-50.0*(alpha/MY_PI - 1.0/2.0)) );
} else { // the central particle is i
if (j < atom->nlocal || k < atom->nlocal) {
const double enx_ij = delx_ij * rinv_ij;
const double eny_ij = dely_ij * rinv_ij;
const double enz_ij = delz_ij * rinv_ij;
const double enx_ik = delx_ik * rinv_ik;
const double eny_ik = dely_ik * rinv_ik;
const double enz_ik = delz_ik * rinv_ik;
const double alpha = std::acos(enx_ij*enx_ik + eny_ij*eny_ik + enz_ij*enz_ik);
pjk[0] += 1.0/( 1.0 + std::exp(-50.0*(alpha/MY_PI - 1.0/2.0)) );
}
std::vector<double> distances = {r_ij,r_ik,r_jk};
auto maxElement = std::max_element(distances.begin(), distances.end());
double maxValue = *maxElement;
int maxIndex = std::distance(distances.begin(), maxElement);
if (maxIndex == 0) { // the central particle is k
const double enx_ki = -delx_ik * rinv_ik;
const double eny_ki = -dely_ik * rinv_ik;
const double enz_ki = -delz_ik * rinv_ik;
const double enx_kj = -delx_jk * rinv_jk;
const double eny_kj = -dely_jk * rinv_jk;
const double enz_kj = -delz_jk * rinv_jk;
const double alpha = std::acos(enx_ki*enx_kj + eny_ki*eny_kj + enz_ki*enz_kj);
pij[0] += 1.0/( 1.0 + std::exp(-50.0*(alpha/MY_PI - 1.0/2.0)));
} else if (maxIndex == 1) { // the central particle is j
const double enx_ji = -delx_ij * rinv_ij;
const double eny_ji = -dely_ij * rinv_ij;
const double enz_ji = -delz_ij * rinv_ij;
const double enx_jk = delx_jk * rinv_jk;
const double eny_jk = dely_jk * rinv_jk;
const double enz_jk = delz_jk * rinv_jk;
const double alpha = std::acos(enx_ji*enx_jk + eny_ji*eny_jk + enz_ji*enz_jk);
pik[0] += 1.0/( 1.0 + std::exp(-50.0*(alpha/MY_PI - 1.0/2.0)));
} else { // the central particle is i
if (j < atom->nlocal || k < atom->nlocal) {
const double enx_ij = delx_ij * rinv_ij;
const double eny_ij = dely_ij * rinv_ij;
const double enz_ij = delz_ij * rinv_ij;
const double enx_ik = delx_ik * rinv_ik;
const double eny_ik = dely_ik * rinv_ik;
const double enz_ik = delz_ik * rinv_ik;
const double alpha = std::acos(enx_ij*enx_ik + eny_ij*eny_ik + enz_ij*enz_ik);
pjk[0] += 1.0/( 1.0 + std::exp(-50.0*(alpha/MY_PI - 1.0/2.0)));
}
}
}
}
}
}
}
}
}
@ -511,16 +486,12 @@ void FixGranularMDR::mean_surf_disp()
NeighList * list = pair->list;
const int size_history = pair->get_size_history();
int i,j,k,ii,jj,inum,jnum,itype,jtype;
int *ilist,*jlist,*numneigh,**firstneigh;
int *touch,**firsttouch;
double *history,*allhistory,**firsthistory;
bool touchflag = false;
class GranularModel* model;
class GranularModel** models_list = pair->models_list;
int ** types_indices = pair->types_indices;
@ -530,7 +501,6 @@ void FixGranularMDR::mean_surf_disp()
double *radius = atom->radius;
int nlocal = atom->nlocal;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
@ -657,11 +627,9 @@ void FixGranularMDR::update_fix_gran_wall(Fix* fix_in)
Region * region = fix->region;
const int size_history = model->size_history;
int i, m, nc, iwall;
double vwall[3];
bool touchflag = false;
int history_update = 1;
model->history_update = history_update;
@ -688,21 +656,21 @@ void FixGranularMDR::update_fix_gran_wall(Fix* fix_in)
nc = region->surface(x[i][0], x[i][1], x[i][2], radius[i] + model->pulloff_distance(radius[i], 0.0));
if (nc == 0) {
fix->ncontact[i] = 0;
continue;
}
if (nc == 1) {
fix->c2r[0] = 0;
iwall = region->contact[0].iwall;
if (fix->ncontact[i] == 0) {
fix->ncontact[i] = 1;
fix->walls[i][0] = iwall;
for (m = 0; m < size_history; m++) fix->history_many[i][0][m] = 0.0;
} else if (fix->ncontact[i] > 1 || iwall != fix->walls[i][0])
fix->update_contacts(i, nc);
} else
if (nc == 0) {
fix->ncontact[i] = 0;
continue;
}
if (nc == 1) {
fix->c2r[0] = 0;
iwall = region->contact[0].iwall;
if (fix->ncontact[i] == 0) {
fix->ncontact[i] = 1;
fix->walls[i][0] = iwall;
for (m = 0; m < size_history; m++) fix->history_many[i][0][m] = 0.0;
} else if (fix->ncontact[i] > 1 || iwall != fix->walls[i][0])
fix->update_contacts(i, nc);
} else
fix->update_contacts(i, nc);
// process current contacts

46
src/GRANULAR/fix_granular_mdr.h
View File

@ -26,44 +26,46 @@ namespace LAMMPS_NS {
class FixGranularMDR : public Fix {
public:
double * Ro;
double * Vgeo;
double * Velas;
double * Vcaps;
double * eps_bar;
double * dRnumerator;
double * dRdenominator;
double * Acon0;
double * Acon1;
double * Atot;
double * Atot_sum;
double * ddelta_bar;
double * psi;
double * psi_b;
double * sigmaxx;
double * sigmayy;
double * sigmazz;
double * history_setup_flag;
double * contacts;
double * adhesive_length;
FixGranularMDR(class LAMMPS *, int, char **);
~FixGranularMDR() override;
int setmask() override;
void post_constructor() override;
void setup(int) override;
void setup_pre_force(int) override;
void pre_force(int) override;
void end_of_step() override; // FOR MDR
void end_of_step() override;
int pack_forward_comm(int, int *, double *, int, int *) override;
void unpack_forward_comm(int, int, double *) override;
private:
int comm_stage;
char *id_fix;
void radius_update();
void mean_surf_disp();
void calculate_contact_penalty();
void update_fix_gran_wall(Fix*);
double *Ro; // initial radius
double *Vgeo; // geometric particle volume of apparent particle afterremoving spherical cap volume
double *Velas; // particle volume from linear elasticity
double *Vcaps; // spherical cap volume from intersection of apparentradius particle and contact planes
double *eps_bar; // volume-averaged infinitesimal strain tensor
double *dRnumerator; // summation of numerator terms in calculation of dR
double *dRdenominator; // summation of denominator terms in calculation of dR
double *Acon0; // total area involved in contacts: Acon^{n}
double *Acon1; // total area involved in contacts: Acon^{n+1}
double *Atot; // total particle area
double *Atot_sum; // running sum of contact area minus cap area
double *ddelta_bar; // change in mean surface displacement
double *psi; // ratio of free surface area to total surface area
double *psi_b; // TEMPORARY, SINCE PSI_B IS ALREADY DEFINED IN THEINPUT SCRIPT
double *sigmaxx; // xx-component of the stress tensor, not necessary forforce calculation
double *sigmayy; // yy-component of the stress tensor, not necessary forforce calculation
double *sigmazz; // zz-component of the stress tensor, not necessary forforce calculation
double *history_setup_flag; // flag to check if history variables have beeninitialized
double *contacts; // total contacts on particle
double *adhesive_length; // total length of adhesive contact on a particle
};
} // namespace LAMMPS_NS

105
src/GRANULAR/gran_sub_mod_normal.cpp
View File

@ -13,12 +13,13 @@
#include "gran_sub_mod_normal.h"
#include "atom.h"
#include "error.h"
#include "citeme.h"
#include "granular_model.h"
#include "math_const.h"
#include "atom.h"
#include "modify.h"
#include "update.h"
#include "citeme.h"
#include <cmath>
#include <iomanip>
@ -414,6 +415,7 @@ GranSubModNormalMDR::GranSubModNormalMDR(GranularModel *gm, LAMMPS *lmp) :
num_coeffs = 6; // Young's Modulus, Poisson's ratio, yield stress, effective surface energy, psi_b, coefficent of restitution
contact_radius_flag = 1;
size_history = 26;
fix_mdr_flag = 0;
nondefault_history_transfer = 1;
transfer_history_factor = new double[size_history];
@ -425,6 +427,14 @@ GranSubModNormalMDR::GranSubModNormalMDR(GranularModel *gm, LAMMPS *lmp) :
/* ---------------------------------------------------------------------- */
GranSubModNormalMDR::~GranSubModNormalMDR()
{
if (fix_mdr_flag)
modify->delete_fix("MDR");
}
/* ---------------------------------------------------------------------- */
void GranSubModNormalMDR::coeffs_to_local()
{
E = coeffs[0]; // Young's modulus
@ -444,6 +454,57 @@ void GranSubModNormalMDR::coeffs_to_local()
/* ---------------------------------------------------------------------- */
void GranSubModNormalMDR::init()
{
if (!fix_mdr_flag) {
if (modify->get_fix_by_style("MDR").size() == 0)
modify->add_fix("MDR all GRANULAR/MDR");
fix_mdr_flag = 1;
}
// initialize particle history variables
int tmp1, tmp2;
int index_Ro = atom->find_custom("Ro", tmp1, tmp2); // initial radius
int index_Vcaps = atom->find_custom("Vcaps", tmp1, tmp2); // spherical cap volume from intersection of apparent radius particle and contact planes
int index_Vgeo = atom->find_custom("Vgeo", tmp1, tmp2); // geometric particle volume of apparent particle after removing spherical cap volume
int index_Velas = atom->find_custom("Velas", tmp1, tmp2); // particle volume from linear elasticity
int index_eps_bar = atom->find_custom("eps_bar", tmp1, tmp2); // volume-averaged infinitesimal strain tensor
int index_dRnumerator = atom->find_custom("dRnumerator", tmp1, tmp2); // summation of numerator terms in calculation of dR
int index_dRdenominator = atom->find_custom("dRdenominator", tmp1, tmp2); // summation of denominator terms in calculation of dR
int index_Acon0 = atom->find_custom("Acon0", tmp1, tmp2); // total area involved in contacts: Acon^{n}
int index_Acon1 = atom->find_custom("Acon1", tmp1, tmp2); // total area involved in contacts: Acon^{n+1}
int index_Atot = atom->find_custom("Atot", tmp1, tmp2); // total particle area
int index_Atot_sum = atom->find_custom("Atot_sum", tmp1, tmp2); // running sum of contact area minus cap area
int index_ddelta_bar = atom->find_custom("ddelta_bar", tmp1, tmp2); // change in mean surface displacement
int index_psi = atom->find_custom("psi", tmp1, tmp2); // ratio of free surface area to total surface area
int index_sigmaxx = atom->find_custom("sigmaxx", tmp1, tmp2); // xx-component of the stress tensor, not necessary for force calculation
int index_sigmayy = atom->find_custom("sigmayy", tmp1, tmp2); // yy-component of the stress tensor, not necessary for force calculation
int index_sigmazz = atom->find_custom("sigmazz", tmp1, tmp2); // zz-component of the stress tensor, not necessary for force calculation
int index_contacts = atom->find_custom("contacts", tmp1, tmp2); // total contacts on particle
int index_adhesive_length = atom->find_custom("adhesive_length", tmp1, tmp2); // total contacts on particle
Rinitial = atom->dvector[index_Ro];
Vgeo = atom->dvector[index_Vgeo];
Velas = atom->dvector[index_Velas];
Vcaps = atom->dvector[index_Vcaps];
eps_bar = atom->dvector[index_eps_bar];
dRnumerator = atom->dvector[index_dRnumerator];
dRdenominator = atom->dvector[index_dRdenominator];
Acon0 = atom->dvector[index_Acon0];
Acon1 = atom->dvector[index_Acon1];
Atot = atom->dvector[index_Atot];
Atot_sum = atom->dvector[index_Atot_sum];
ddelta_bar = atom->dvector[index_ddelta_bar];
psi = atom->dvector[index_psi];
sigmaxx = atom->dvector[index_sigmaxx];
sigmayy = atom->dvector[index_sigmayy];
sigmazz = atom->dvector[index_sigmazz];
contacts = atom->dvector[index_contacts];
adhesive_length = atom->dvector[index_adhesive_length];
}
/* ---------------------------------------------------------------------- */
double GranSubModNormalMDR::calculate_forces()
{
@ -507,46 +568,6 @@ double GranSubModNormalMDR::calculate_forces()
const int deltap_offset_0 = 24;
const int deltap_offset_1 = 25;
// initialize particle history variables
int tmp1, tmp2;
int index_Ro = atom->find_custom("Ro",tmp1,tmp2); // initial radius
int index_Vcaps = atom->find_custom("Vcaps",tmp1,tmp2); // spherical cap volume from intersection of apparent radius particle and contact planes
int index_Vgeo = atom->find_custom("Vgeo",tmp1,tmp2); // geometric particle volume of apparent particle after removing spherical cap volume
int index_Velas = atom->find_custom("Velas",tmp1,tmp2); // particle volume from linear elasticity
int index_eps_bar = atom->find_custom("eps_bar",tmp1,tmp2); // volume-averaged infinitesimal strain tensor
int index_dRnumerator = atom->find_custom("dRnumerator",tmp1,tmp2); // summation of numerator terms in calculation of dR
int index_dRdenominator = atom->find_custom("dRdenominator",tmp1,tmp2); // summation of denominator terms in calculation of dR
int index_Acon0 = atom->find_custom("Acon0",tmp1,tmp2); // total area involved in contacts: Acon^{n}
int index_Acon1 = atom->find_custom("Acon1",tmp1,tmp2); // total area involved in contacts: Acon^{n+1}
int index_Atot = atom->find_custom("Atot",tmp1,tmp2); // total particle area
int index_Atot_sum = atom->find_custom("Atot_sum",tmp1,tmp2); // running sum of contact area minus cap area
int index_ddelta_bar = atom->find_custom("ddelta_bar",tmp1,tmp2); // change in mean surface displacement
int index_psi = atom->find_custom("psi",tmp1,tmp2); // ratio of free surface area to total surface area
int index_sigmaxx = atom->find_custom("sigmaxx",tmp1,tmp2); // xx-component of the stress tensor, not necessary for force calculation
int index_sigmayy = atom->find_custom("sigmayy",tmp1,tmp2); // yy-component of the stress tensor, not necessary for force calculation
int index_sigmazz = atom->find_custom("sigmazz",tmp1,tmp2); // zz-component of the stress tensor, not necessary for force calculation
int index_contacts = atom->find_custom("contacts",tmp1,tmp2); // total contacts on particle
int index_adhesive_length = atom->find_custom("adhesive_length",tmp1,tmp2); // total contacts on particle
double * Rinitial = atom->dvector[index_Ro];
double * Vgeo = atom->dvector[index_Vgeo];
double * Velas = atom->dvector[index_Velas];
double * Vcaps = atom->dvector[index_Vcaps];
double * eps_bar = atom->dvector[index_eps_bar];
double * dRnumerator = atom->dvector[index_dRnumerator];
double * dRdenominator = atom->dvector[index_dRdenominator];
double * Acon0 = atom->dvector[index_Acon0];
double * Acon1 = atom->dvector[index_Acon1];
double * Atot = atom->dvector[index_Atot];
double * Atot_sum = atom->dvector[index_Atot_sum];
double * ddelta_bar = atom->dvector[index_ddelta_bar];
double * psi = atom->dvector[index_psi];
double * sigmaxx = atom->dvector[index_sigmaxx];
double * sigmayy = atom->dvector[index_sigmayy];
double * sigmazz = atom->dvector[index_sigmazz];
double * contacts = atom->dvector[index_contacts];
double * adhesive_length = atom->dvector[index_adhesive_length];
double * history = & gm->history[history_index]; // load in all history variables
// Rigid flat placement scheme

7
src/GRANULAR/gran_sub_mod_normal.h
View File

@ -139,13 +139,20 @@ namespace Granular_NS {
class GranSubModNormalMDR : public GranSubModNormal {
public:
GranSubModNormalMDR(class GranularModel *, class LAMMPS *);
~GranSubModNormalMDR() override;
void coeffs_to_local() override;
void init() override;
double calculate_forces() override;
void set_fncrit() override;
double psi_b;
protected:
double E, nu, Y, gamma, CoR, F;
double *Rinitial, *Vgeo, *Velas, *Vcaps, *eps_bar, *dRnumerator;
double *dRdenominator, *Acon0, *Acon1, *Atot, *Atot_sum, *ddelta_bar;
double *psi, *sigmaxx, *sigmayy, *sigmazz, *contacts, *adhesive_length;
int fix_mdr_flag;
};
} // namespace Granular_NS

37
src/GRANULAR/pair_granular.cpp
View File

@ -78,8 +78,6 @@ PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp)
fix_history = nullptr;
fix_dummy = dynamic_cast<FixDummy *>(modify->add_fix("NEIGH_HISTORY_GRANULAR_DUMMY all DUMMY"));
fix_flag = 0;
}
/* ---------------------------------------------------------------------- */
@ -475,41 +473,6 @@ void PairGranular::init_style()
if (!fix_history) error->all(FLERR,"Could not find pair fix neigh history ID");
}
if (model->sub_models[NORMAL]->name == "mdr") {
//Store persistent per atom quantities
if (! fix_flag) {
int tmp1, tmp2;
const char * id_fix = "MDR_PARTICLE_HISTORY_VARIABLES";
modify->add_fix(fmt::format("{} all property/atom d_Ro d_Vcaps d_Vgeo d_Velas d_eps_bar d_dRnumerator d_dRdenominator d_Acon0 d_Acon1 d_Atot d_Atot_sum d_ddelta_bar d_psi d_psi_b d_history_setup_flag d_sigmaxx d_sigmayy d_sigmazz d_contacts d_adhesive_length ghost yes", id_fix));
// d2_volSums 4 --> allows an array of 4 to defined.
index_Ro = atom->find_custom("Ro",tmp1,tmp2); // initial radius
index_Vcaps = atom->find_custom("Vcaps",tmp1,tmp2); // spherical cap volume from intersection of apparent radius particle and contact planes
index_Vgeo = atom->find_custom("Vgeo",tmp1,tmp2); // geometric particle volume of apparent particle after removing spherical cap volume
index_Velas = atom->find_custom("Velas",tmp1,tmp2); // particle volume from linear elasticity
index_eps_bar = atom->find_custom("eps_bar",tmp1,tmp2); // volume-averaged infinitesimal strain tensor
index_dRnumerator = atom->find_custom("dRnumerator",tmp1,tmp2); // summation of numerator terms in calculation of dR
index_dRdenominator = atom->find_custom("dRdenominator",tmp1,tmp2); // summation of denominator terms in calculation of dR
index_Acon0 = atom->find_custom("Acon0",tmp1,tmp2); // total area involved in contacts: Acon^{n}
index_Acon1 = atom->find_custom("Acon1",tmp1,tmp2); // total area involved in contacts: Acon^{n+1}
index_Atot = atom->find_custom("Atot",tmp1,tmp2); // total particle area
index_Atot_sum = atom->find_custom("Atot_sum",tmp1,tmp2); // running sum of contact area minus cap area
index_ddelta_bar = atom->find_custom("ddelta_bar",tmp1,tmp2); // change in mean surface displacement
index_psi = atom->find_custom("psi",tmp1,tmp2); // ratio of free surface area to total surface area
index_psi_b = atom->find_custom("psi_b",tmp1,tmp2); // TEMPORARY, SINCE PSI_B IS ALREADY DEFINED IN THE INPUT SCRIPT
index_history_setup_flag = atom->find_custom("history_setup_flag",tmp1,tmp2); // flag to check if history variables have been initialized
index_sigmaxx = atom->find_custom("sigmaxx",tmp1,tmp2); // xx-component of the stress tensor, not necessary for force calculation
index_sigmayy = atom->find_custom("sigmayy",tmp1,tmp2); // yy-component of the stress tensor, not necessary for force calculation
index_sigmazz = atom->find_custom("sigmazz",tmp1,tmp2); // zz-component of the stress tensor, not necessary for force calculation
index_contacts = atom->find_custom("contacts",tmp1,tmp2); // total contacts on particle
index_adhesive_length = atom->find_custom("adhesive_length",tmp1,tmp2); // total length of adhesive contact on a particle
// Initiate MDR radius update fix
modify->add_fix("fix_mdr all GRANULAR/MDR");
fix_flag = 1;
}
}
// check for FixFreeze and set freeze_group_bit
auto fixlist = modify->get_fix_by_style("^freeze");

25
src/GRANULAR/pair_granular.h
View File

@ -53,8 +53,6 @@ class PairGranular : public Pair {
int **types_indices;
int nmodels, maxmodels;
class FixStoreAtom * fix_store;
protected:
int freeze_group_bit;
int use_history;
@ -67,29 +65,6 @@ class PairGranular : public Pair {
class FixDummy *fix_dummy;
class FixNeighHistory *fix_history;
// MDR particle history variables
int fix_flag;
int index_Ro;
int index_Vcaps;
int index_Vgeo;
int index_Velas;
int index_eps_bar;
int index_dRnumerator;
int index_dRdenominator;
int index_Acon0;
int index_Acon1;
int index_Atot;
int index_Atot_sum;
int index_ddelta_bar;
int index_psi;
int index_psi_b;
int index_history_setup_flag;
int index_sigmaxx;
int index_sigmayy;
int index_sigmazz;
int index_contacts;
int index_adhesive_length;
// storage of rigid body masses for use in granular interactions
class Fix *fix_rigid; // ptr to rigid body fix, null pointer if none