From 9a6dc2ff11dd846f4e2cf7b4cada63af21dcb110 Mon Sep 17 00:00:00 2001 From: "Dan S. Bolintineanu" Date: Wed, 6 Mar 2019 13:54:32 -0700 Subject: [PATCH] Removed several files that should not have been included --- src/GRANULAR/pair_gran_dmt_rolling.cpp | 721 -------- src/GRANULAR/pair_gran_dmt_rolling.h | 55 - src/GRANULAR/pair_gran_dmt_rolling2.cpp | 719 -------- .../pair_gran_hooke_history_multi.cpp | 915 ---------- src/GRANULAR/pair_gran_hooke_history_multi.h | 109 -- src/GRANULAR/pair_gran_jkr_rolling.cpp | 763 -------- src/GRANULAR/pair_gran_jkr_rolling.h | 56 - src/GRANULAR/pair_gran_jkr_rolling_multi.cpp | 1181 ------------ src/GRANULAR/pair_gran_jkr_rolling_multi.h | 87 - src/GRANULAR/pair_granular_multi.cpp | 1624 ----------------- src/GRANULAR/pair_granular_multi.h | 107 -- 11 files changed, 6337 deletions(-) delete mode 100644 src/GRANULAR/pair_gran_dmt_rolling.cpp delete mode 100644 src/GRANULAR/pair_gran_dmt_rolling.h delete mode 100644 src/GRANULAR/pair_gran_dmt_rolling2.cpp delete mode 100644 src/GRANULAR/pair_gran_hooke_history_multi.cpp delete mode 100644 src/GRANULAR/pair_gran_hooke_history_multi.h delete mode 100644 src/GRANULAR/pair_gran_jkr_rolling.cpp delete mode 100644 src/GRANULAR/pair_gran_jkr_rolling.h delete mode 100644 src/GRANULAR/pair_gran_jkr_rolling_multi.cpp delete mode 100644 src/GRANULAR/pair_gran_jkr_rolling_multi.h delete mode 100644 src/GRANULAR/pair_granular_multi.cpp delete mode 100644 src/GRANULAR/pair_granular_multi.h diff --git a/src/GRANULAR/pair_gran_dmt_rolling.cpp b/src/GRANULAR/pair_gran_dmt_rolling.cpp deleted file mode 100644 index f293998e92..0000000000 --- a/src/GRANULAR/pair_gran_dmt_rolling.cpp +++ /dev/null @@ -1,721 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. - ------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - Contributing authors: Leo Silbert (SNL), Gary Grest (SNL) - ------------------------------------------------------------------------- */ - -#include -#include -#include -#include -#include "pair_gran_dmt_rolling.h" -#include "atom.h" -#include "update.h" -#include "force.h" -#include "fix.h" -#include "fix_neigh_history.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "comm.h" -#include "memory.h" -#include "error.h" -#include "math_const.h" - -using namespace LAMMPS_NS; -using namespace MathConst; - -#define TWOTHIRDS 0.6666666666666666 -#define EPSILON 1e-10 - -enum {TSUJI, BRILLIANTOV}; -enum {INDEP, BRILLROLL}; - -/* ---------------------------------------------------------------------- */ - -PairGranDMTRolling::PairGranDMTRolling(LAMMPS *lmp) : - PairGranHookeHistory(lmp, 7), - E_one(0), G_one(0), pois(0), muS_one(0), cor(0), alpha_one(0), - Ecoh_one(0), kR_one(0), muR_one(0), etaR_one(0) -{ - int ntypes = atom->ntypes; - memory->create(E,ntypes+1,ntypes+1,"pair:E"); - memory->create(G,ntypes+1,ntypes+1,"pair:G"); - memory->create(alpha,ntypes+1,ntypes+1,"pair:alpha"); - memory->create(gamman,ntypes+1,ntypes+1,"pair:gamman"); - memory->create(muS,ntypes+1,ntypes+1,"pair:muS"); - memory->create(Ecoh,ntypes+1,ntypes+1,"pair:Ecoh"); - memory->create(kR,ntypes+1,ntypes+1,"pair:kR"); - memory->create(muR,ntypes+1,ntypes+1,"pair:muR"); - memory->create(etaR,ntypes+1,ntypes+1,"pair:etaR"); -} - -/* ---------------------------------------------------------------------- */ -PairGranDMTRolling::~PairGranDMTRolling() -{ - delete [] E_one; - delete [] G_one; - delete [] pois; - delete [] muS_one; - delete [] cor; - delete [] alpha_one; - delete [] Ecoh_one; - delete [] kR_one; - delete [] muR_one; - delete [] etaR_one; - //TODO: Make all this work with standard pair coeff type commands. - //Also these should not be in the destructor. - memory->destroy(E); - memory->destroy(G); - memory->destroy(alpha); - memory->destroy(gamman); - memory->destroy(muS); - memory->destroy(Ecoh); - memory->destroy(kR); - memory->destroy(muR); - memory->destroy(etaR); -} -/* ---------------------------------------------------------------------- */ - -void PairGranDMTRolling::compute(int eflag, int vflag) -{ - int i,j,ii,jj,inum,jnum; - int itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,nx,ny,nz; - double radi,radj,radsum,rsq,r,rinv,rsqinv,R,a; - double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; - double wr1,wr2,wr3; - double vtr1,vtr2,vtr3,vrel; - double kn, kt, k_Q, k_R, eta_N, eta_T, eta_Q, eta_R; - double Fhz, Fdamp, Fdmt, Fne, Fntot, Fscrit, Frcrit; - double overlap; - double mi,mj,meff,damp,ccel,tor1,tor2,tor3; - double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3,vrlmag,vrlmaginv; - double rollmag, rolldotn, scalefac; - double fr, fr1, fr2, fr3; - double signtwist, magtwist, magtortwist, Mtcrit; - double fs,fs1,fs2,fs3,roll1,roll2,roll3,torroll1,torroll2,torroll3; - double tortwist1, tortwist2, tortwist3; - double shrmag,rsht; - int *ilist,*jlist,*numneigh,**firstneigh; - int *touch,**firsttouch; - double *shear,*allshear,**firstshear; - - if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = vflag_fdotr = 0; - - int shearupdate = 1; - if (update->setupflag) shearupdate = 0; - - // update rigid body info for owned & ghost atoms if using FixRigid masses - // body[i] = which body atom I is in, -1 if none - // mass_body = mass of each rigid body - - if (fix_rigid && neighbor->ago == 0){ - int tmp; - int *body = (int *) fix_rigid->extract("body",tmp); - double *mass_body = (double *) fix_rigid->extract("masstotal",tmp); - if (atom->nmax > nmax) { - memory->destroy(mass_rigid); - nmax = atom->nmax; - memory->create(mass_rigid,nmax,"pair:mass_rigid"); - } - int nlocal = atom->nlocal; - for (i = 0; i < nlocal; i++) - if (body[i] >= 0) mass_rigid[i] = mass_body[body[i]]; - else mass_rigid[i] = 0.0; - comm->forward_comm_pair(this); - } - - double **x = atom->x; - double **v = atom->v; - double **f = atom->f; - double **omega = atom->omega; - double **torque = atom->torque; - double *radius = atom->radius; - double *rmass = atom->rmass; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - inum = list->inum; - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->firstneigh; - firsttouch = fix_history->firstflag; - firstshear = fix_history->firstvalue; - - // loop over neighbors of my atoms - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - itype = type[i]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - radi = radius[i]; - touch = firsttouch[i]; - allshear = firstshear[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - jtype = type[j]; - j &= NEIGHMASK; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - radj = radius[j]; - radsum = radi + radj; - - if (rsq >= radsum*radsum){ - // unset non-touching neighbors - touch[jj] = 0; - shear = &allshear[size_history*jj]; - for (int k = 0; k < size_history; k++) - shear[k] = 0.0; - } else { - r = sqrt(rsq); - rinv = 1.0/r; - rsqinv = 1.0/rsq; - R = radi*radj/(radi+radj); - nx = delx*rinv; - ny = dely*rinv; - nz = delz*rinv; - - // relative translational velocity - - vr1 = v[i][0] - v[j][0]; - vr2 = v[i][1] - v[j][1]; - vr3 = v[i][2] - v[j][2]; - - // normal component - - vnnr = vr1*nx + vr2*ny + vr3*nz; //v_R . n - vn1 = nx*vnnr; - vn2 = ny*vnnr; - vn3 = nz*vnnr; - - // meff = effective mass of pair of particles - // if I or J part of rigid body, use body mass - // if I or J is frozen, meff is other particle - - mi = rmass[i]; - mj = rmass[j]; - if (fix_rigid) { - if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; - if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; - } - - meff = mi*mj / (mi+mj); - if (mask[i] & freeze_group_bit) meff = mj; - if (mask[j] & freeze_group_bit) meff = mi; - - //**************************************** - //Normal force = Hertzian contact + DMT + damping - //**************************************** - overlap = radsum - r; - a = sqrt(R*overlap); - kn = 4.0/3.0*E[itype][jtype]*a; - Fhz = kn*overlap; - - //Damping (based on Tsuji et al) - if (normaldamp == BRILLIANTOV) eta_N = a*meff*gamman[itype][jtype]; - else if (normaldamp == TSUJI) eta_N=alpha[itype][jtype]*sqrt(meff*kn); - - Fdamp = -eta_N*vnnr; //F_nd eq 23 and Zhao eq 19 - - //DMT - Fdmt = -4*MY_PI*Ecoh[itype][jtype]*R; - - Fne = Fhz + Fdmt; - Fntot = Fne + Fdamp; - - //**************************************** - //Tangential force, including shear history effects - //**************************************** - - // tangential component - vt1 = vr1 - vn1; - vt2 = vr2 - vn2; - vt3 = vr3 - vn3; - - // relative rotational velocity - //Luding Gran Matt 2008, v10,p235 suggests correcting radi and radj by subtracting - //delta/2, i.e. instead of radi, use distance to center of contact point? - wr1 = (radi*omega[i][0] + radj*omega[j][0]); - wr2 = (radi*omega[i][1] + radj*omega[j][1]); - wr3 = (radi*omega[i][2] + radj*omega[j][2]); - - // relative tangential velocities - vtr1 = vt1 - (nz*wr2-ny*wr3); - vtr2 = vt2 - (nx*wr3-nz*wr1); - vtr3 = vt3 - (ny*wr1-nx*wr2); - vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; - vrel = sqrt(vrel); - - // shear history effects - touch[jj] = 1; - shear = &allshear[size_history*jj]; - shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + - shear[2]*shear[2]); - - // Rotate and update shear displacements. - // See e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 - if (shearupdate) { - rsht = shear[0]*nx + shear[1]*ny + shear[2]*nz; - if (fabs(rsht) < EPSILON) rsht = 0; - if (rsht > 0){ - scalefac = shrmag/(shrmag - rsht); //if rhst == shrmag, contacting pair has rotated 90 deg. in one step, in which case you deserve a crash! - shear[0] -= rsht*nx; - shear[1] -= rsht*ny; - shear[2] -= rsht*nz; - //Also rescale to preserve magnitude - shear[0] *= scalefac; - shear[1] *= scalefac; - shear[2] *= scalefac; - } - //Update shear history - shear[0] += vtr1*dt; - shear[1] += vtr2*dt; - shear[2] += vtr3*dt; - } - - // tangential forces = shear + tangential velocity damping - // following Zhao and Marshall Phys Fluids v20, p043302 (2008) - kt=8.0*G[itype][jtype]*a; - - eta_T = eta_N; //Based on discussion in Marshall; eta_T can also be an independent parameter - fs1 = -kt*shear[0] - eta_T*vtr1; //eq 26 - fs2 = -kt*shear[1] - eta_T*vtr2; - fs3 = -kt*shear[2] - eta_T*vtr3; - - // rescale frictional displacements and forces if needed - Fscrit = muS[itype][jtype] * fabs(Fne); - // For JKR, use eq 43 of Marshall. For DMT, use Fne instead - shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + - shear[2]*shear[2]); - fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); - if (fs > Fscrit) { - if (shrmag != 0.0) { - //shear[0] = (Fcrit/fs) * (shear[0] + eta_T*vtr1/kt) - eta_T*vtr1/kt; - //shear[1] = (Fcrit/fs) * (shear[1] + eta_T*vtr1/kt) - eta_T*vtr1/kt; - //shear[2] = (Fcrit/fs) * (shear[2] + eta_T*vtr1/kt) - eta_T*vtr1/kt; - shear[0] = -1.0/kt*(Fscrit*fs1/fs + eta_T*vtr1); //Same as above, but simpler (check!) - shear[1] = -1.0/kt*(Fscrit*fs2/fs + eta_T*vtr2); - shear[2] = -1.0/kt*(Fscrit*fs3/fs + eta_T*vtr3); - fs1 *= Fscrit/fs; - fs2 *= Fscrit/fs; - fs3 *= Fscrit/fs; - } else fs1 = fs2 = fs3 = 0.0; - } - - //**************************************** - // Rolling force, including shear history effects - //**************************************** - - relrot1 = omega[i][0] - omega[j][0]; - relrot2 = omega[i][1] - omega[j][1]; - relrot3 = omega[i][2] - omega[j][2]; - - // rolling velocity, see eq. 31 of Wang et al, Particuology v 23, p 49 (2015) - // This is different from the Marshall papers, which use the Bagi/Kuhn formulation - // for rolling velocity (see Wang et al for why the latter is wrong) - vrl1 = R*(relrot2*nz - relrot3*ny); //- 0.5*((radj-radi)/radsum)*vtr1; - vrl2 = R*(relrot3*nx - relrot1*nz); //- 0.5*((radj-radi)/radsum)*vtr2; - vrl3 = R*(relrot1*ny - relrot2*nx); //- 0.5*((radj-radi)/radsum)*vtr3; - vrlmag = sqrt(vrl1*vrl1+vrl2*vrl2+vrl3*vrl3); - if (vrlmag != 0.0) vrlmaginv = 1.0/vrlmag; - else vrlmaginv = 0.0; - - // Rolling displacement - rollmag = sqrt(shear[3]*shear[3] + shear[4]*shear[4] + shear[5]*shear[5]); - rolldotn = shear[3]*nx + shear[4]*ny + shear[5]*nz; - - if (shearupdate) { - if (fabs(rolldotn) < EPSILON) rolldotn = 0; - if (rolldotn > 0){ //Rotate into tangential plane - scalefac = rollmag/(rollmag - rolldotn); - shear[3] -= rolldotn*nx; - shear[4] -= rolldotn*ny; - shear[5] -= rolldotn*nz; - //Also rescale to preserve magnitude - shear[3] *= scalefac; - shear[4] *= scalefac; - shear[5] *= scalefac; - } - shear[3] += vrl1*dt; - shear[4] += vrl2*dt; - shear[5] += vrl3*dt; - } - - k_R = kR[itype][jtype]; - if (rollingdamp == INDEP) eta_R = etaR[itype][jtype]; - else if (rollingdamp == BRILLROLL) eta_R = muR[itype][jtype]*fabs(Fne); - fr1 = -k_R*shear[3] - eta_R*vrl1; - fr2 = -k_R*shear[4] - eta_R*vrl2; - fr3 = -k_R*shear[5] - eta_R*vrl3; - - // rescale frictional displacements and forces if needed - Frcrit = muR[itype][jtype] * fabs(Fne); - - rollmag = sqrt(shear[3]*shear[3] + shear[4]*shear[4] + shear[5]*shear[5]); - fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3); - if (fr > Frcrit) { - if (rollmag != 0.0) { - shear[3] = -1.0/k_R*(Frcrit*fr1/fr + eta_R*vrl1); - shear[4] = -1.0/k_R*(Frcrit*fr2/fr + eta_R*vrl2); - shear[5] = -1.0/k_R*(Frcrit*fr3/fr + eta_R*vrl3); - fr1 *= Frcrit/fr; - fr2 *= Frcrit/fr; - fr3 *= Frcrit/fr; - } else fr1 = fr2 = fr3 = 0.0; - } - - - //**************************************** - // Twisting torque, including shear history effects - //**************************************** - magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of Marshall) - shear[6] += magtwist*dt; - k_Q = 0.5*kt*a*a;; //eq 32 - eta_Q = 0.5*eta_T*a*a; - magtortwist = -k_Q*shear[6] - eta_Q*magtwist;//M_t torque (eq 30) - - signtwist = (magtwist > 0) - (magtwist < 0); - Mtcrit=TWOTHIRDS*a*Fscrit;//critical torque (eq 44) - if (fabs(magtortwist) > Mtcrit){ - shear[6] = 1.0/k_Q*(Mtcrit*signtwist - eta_Q*magtwist); - magtortwist = -Mtcrit * signtwist; //eq 34 - } - - // Apply forces & torques - - fx = nx*Fntot + fs1; - fy = ny*Fntot + fs2; - fz = nz*Fntot + fs3; - - f[i][0] += fx; - f[i][1] += fy; - f[i][2] += fz; - - tor1 = ny*fs3 - nz*fs2; - tor2 = nz*fs1 - nx*fs3; - tor3 = nx*fs2 - ny*fs1; - - torque[i][0] -= radi*tor1; - torque[i][1] -= radi*tor2; - torque[i][2] -= radi*tor3; - - tortwist1 = magtortwist * nx; - tortwist2 = magtortwist * ny; - tortwist3 = magtortwist * nz; - - torque[i][0] += tortwist1; - torque[i][1] += tortwist2; - torque[i][2] += tortwist3; - - torroll1 = R*(ny*fr3 - nz*fr2); //n cross fr - torroll2 = R*(nz*fr1 - nx*fr3); - torroll3 = R*(nx*fr2 - ny*fr1); - - torque[i][0] += torroll1; - torque[i][1] += torroll2; - torque[i][2] += torroll3; - - if (force->newton_pair || j < nlocal) { - f[j][0] -= fx; - f[j][1] -= fy; - f[j][2] -= fz; - - torque[j][0] -= radj*tor1; - torque[j][1] -= radj*tor2; - torque[j][2] -= radj*tor3; - - torque[j][0] -= tortwist1; - torque[j][1] -= tortwist2; - torque[j][2] -= tortwist3; - - torque[j][0] -= torroll1; - torque[j][1] -= torroll2; - torque[j][2] -= torroll3; - } - if (evflag) ev_tally_xyz(i,j,nlocal,0, - 0.0,0.0,fx,fy,fz,delx,dely,delz); - } - } - } -} - -/* ---------------------------------------------------------------------- - global settings - ------------------------------------------------------------------------- */ - -void PairGranDMTRolling::settings(int narg, char **arg) -{ - if (narg < 6) error->all(FLERR,"Illegal pair_style command"); - - int ntypes = atom->ntypes; - - if (narg < 8*ntypes) error->all(FLERR,"Illegal pair_style command"); - - E_one = new double[ntypes+1]; - G_one = new double[ntypes+1]; - pois = new double[ntypes+1]; - muS_one = new double[ntypes+1]; - cor = new double[ntypes+1]; - alpha_one = new double[ntypes+1]; - Ecoh_one = new double[ntypes+1]; - kR_one = new double[ntypes+1]; - muR_one = new double[ntypes+1]; - etaR_one = new double[ntypes+1]; - - for (int i=0; i < ntypes;i++){ - E_one[i+1] = force->numeric(FLERR, arg[i]); - G_one[i+1] = force->numeric(FLERR, arg[ntypes+i]); - muS_one[i+1] = force->numeric(FLERR, arg[2*ntypes+i]); - cor[i+1] = force->numeric(FLERR, arg[3*ntypes+i]); - Ecoh_one[i+1] = force->numeric(FLERR, arg[4*ntypes+i]); - kR_one[i+1] = force->numeric(FLERR, arg[5*ntypes+i]); - muR_one[i+1] = force->numeric(FLERR, arg[6*ntypes+i]); - etaR_one[i+1] = force->numeric(FLERR, arg[7*ntypes+i]); - } - - //Defaults - normaldamp = TSUJI; - rollingdamp = INDEP; - - int iarg = 8*ntypes; - while (iarg < narg){ - if (strcmp(arg[iarg],"normaldamp") == 0){ - if (iarg+2 > narg) error->all(FLERR, "Invalid pair/gran/dmt/rolling entry"); - if (strcmp(arg[iarg+1],"tsuji") == 0) normaldamp = TSUJI; - else if (strcmp(arg[iarg+1],"brilliantov") == 0) normaldamp = BRILLIANTOV; - else error->all(FLERR, "Invalid normal damping model for pair/gran/dmt/rolling"); - iarg += 2; - } - else if (strcmp(arg[iarg],"rollingdamp") == 0){ - if (iarg+2 > narg) error->all(FLERR, "Invalid pair/gran/dmt/rolling entry"); - if (strcmp(arg[iarg+1],"independent") == 0) rollingdamp = INDEP; - else if (strcmp(arg[iarg+1],"brilliantov") == 0) rollingdamp = BRILLROLL; - else error->all(FLERR, "Invalid rolling damping model for pair/gran/dmt/rolling"); - iarg += 2; - } - else{ - iarg +=1; - } - } - - //Derived from inputs - for (int i=1; i <= ntypes; i++){ - pois[i] = E_one[i]/(2.0*G_one[i]) - 1.0; - alpha_one[i] = 1.2728-4.2783*cor[i]+11.087*cor[i]*cor[i]-22.348*cor[i]*cor[i]*cor[i]+27.467*cor[i]*cor[i]*cor[i]*cor[i]-18.022*cor[i]*cor[i]*cor[i]*cor[i]*cor[i]+4.8218*cor[i]*cor[i]*cor[i]*cor[i]*cor[i]*cor[i]; - for (int j=i; j <= ntypes; j++){ - E[i][j] = E[j][i] = 1/((1-pois[i]*pois[i])/E_one[i]+(1-pois[j]*pois[j])/E_one[j]); - G[i][j] = G[j][i] = 1/((2-pois[i])/G_one[i]+(2-pois[j])/G_one[j]); - if (normaldamp == TSUJI){ - alpha[i][j] = alpha[j][i] = sqrt(alpha_one[i]*alpha_one[j]); - } - else if (normaldamp == BRILLIANTOV){ - gamman[i][j] = gamman[j][i] = sqrt(cor[i]*cor[j]); - } - muS[i][j] = muS[j][i] = sqrt(muS_one[i]*muS_one[j]); - Ecoh[i][j] = Ecoh[j][i] = sqrt(Ecoh_one[i]*Ecoh_one[j]); - kR[i][j] = kR[j][i] = sqrt(kR_one[i]*kR_one[j]); - etaR[i][j] = etaR[j][i] = sqrt(etaR_one[i]*etaR_one[j]); - muR[i][j] = muR[j][i] = sqrt(muR_one[i]*muR_one[j]); - } - } -} - -/* ---------------------------------------------------------------------- */ - -double PairGranDMTRolling::single(int i, int j, int itype, int jtype, - double rsq, - double factor_coul, double factor_lj, - double &fforce) -{ - double radi,radj,radsum; - double r,rinv,rsqinv,delx,dely,delz, nx, ny, nz, R; - double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3,wr1,wr2,wr3; - double overlap, a; - double mi,mj,meff,damp,kn,kt; - double Fhz,Fdamp,Fdmt,Fne,Fntot,Fscrit; - double eta_N,eta_T; - double vtr1,vtr2,vtr3,vrel; - double fs1,fs2,fs3,fs; - double shrmag; - - - double *radius = atom->radius; - radi = radius[i]; - radj = radius[j]; - radsum = radi + radj; - - if (rsq >= radsum*radsum) { - fforce = 0.0; - svector[0] = svector[1] = svector[2] = svector[3] = 0.0; - return 0.0; - } - - r = sqrt(rsq); - rinv = 1.0/r; - rsqinv = 1.0/rsq; - R = radi*radj/radsum; - - // relative translational velocity - - double **v = atom->v; - vr1 = v[i][0] - v[j][0]; - vr2 = v[i][1] - v[j][1]; - vr3 = v[i][2] - v[j][2]; - - // normal component - - double **x = atom->x; - delx = x[i][0] - x[j][0]; - dely = x[i][1] - x[j][1]; - delz = x[i][2] - x[j][2]; - - nx = delx*rinv; - ny = dely*rinv; - nz = delz*rinv; - - - vnnr = vr1*nx + vr2*ny + vr3*nz; - vn1 = nx*vnnr; - vn2 = ny*vnnr; - vn3 = nz*vnnr; - - // tangential component - - vt1 = vr1 - vn1; - vt2 = vr2 - vn2; - vt3 = vr3 - vn3; - - // relative rotational velocity - - double **omega = atom->omega; - wr1 = (radi*omega[i][0] + radj*omega[j][0]); - wr2 = (radi*omega[i][1] + radj*omega[j][1]); - wr3 = (radi*omega[i][2] + radj*omega[j][2]); - - // meff = effective mass of pair of particles - // if I or J part of rigid body, use body mass - // if I or J is frozen, meff is other particle - - double *rmass = atom->rmass; - int *type = atom->type; - int *mask = atom->mask; - - mi = rmass[i]; - mj = rmass[j]; - if (fix_rigid) { - // NOTE: ensure mass_rigid is current for owned+ghost atoms? - if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; - if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; - } - - meff = mi*mj / (mi+mj); - if (mask[i] & freeze_group_bit) meff = mj; - if (mask[j] & freeze_group_bit) meff = mi; - - - // normal force = Hertzian contact + normal velocity damping - overlap = radsum - r; - a = sqrt(R*overlap); - kn = 4.0/3.0*E[itype][jtype]*a; - Fhz = kn*overlap; - - //Damping (based on Tsuji et al) - - eta_N=alpha[itype][jtype]*sqrt(meff*kn); - Fdamp = -eta_N*vnnr; //F_nd eq 23 and Zhao eq 19 - - //DMT - Fdmt = -4*MY_PI*Ecoh[itype][jtype]*R; - - Fne = Fhz + Fdmt; - Fntot = Fne + Fdamp; - - // relative velocities - - vtr1 = vt1 - (nz*wr2-ny*wr3); - vtr2 = vt2 - (nx*wr3-nz*wr1); - vtr3 = vt3 - (ny*wr1-nx*wr2); - vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; - vrel = sqrt(vrel); - - // shear history effects - // neighprev = index of found neigh on previous call - // search entire jnum list of neighbors of I for neighbor J - // start from neighprev, since will typically be next neighbor - // reset neighprev to 0 as necessary - - int jnum = list->numneigh[i]; - int *jlist = list->firstneigh[i]; - double *allshear = fix_history->firstvalue[i]; - - for (int jj = 0; jj < jnum; jj++) { - neighprev++; - if (neighprev >= jnum) neighprev = 0; - if (jlist[neighprev] == j) break; - } - - double *shear = &allshear[size_history*neighprev]; - shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + - shear[2]*shear[2]); - - // tangential forces = shear + tangential velocity damping - kt=8.0*G[itype][jtype]*a; - - eta_T = eta_N; - fs1 = -kt*shear[0] - eta_T*vtr1; - fs2 = -kt*shear[1] - eta_T*vtr2; - fs3 = -kt*shear[2] - eta_T*vtr3; - - // rescale frictional displacements and forces if needed - - fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); - Fscrit= muS[itype][jtype] * fabs(Fne); - - if (fs > Fscrit) { - if (shrmag != 0.0) { - fs1 *= Fscrit/fs; - fs2 *= Fscrit/fs; - fs3 *= Fscrit/fs; - fs *= Fscrit/fs; - } else fs1 = fs2 = fs3 = fs = 0.0; - } - - // set all forces and return no energy - - fforce = Fntot; - - // set single_extra quantities - - svector[0] = fs1; - svector[1] = fs2; - svector[2] = fs3; - svector[3] = fs; - svector[4] = vn1; - svector[5] = vn2; - svector[6] = vn3; - svector[7] = vt1; - svector[8] = vt2; - svector[9] = vt3; - return 0.0; -} diff --git a/src/GRANULAR/pair_gran_dmt_rolling.h b/src/GRANULAR/pair_gran_dmt_rolling.h deleted file mode 100644 index 8f4ae2005e..0000000000 --- a/src/GRANULAR/pair_gran_dmt_rolling.h +++ /dev/null @@ -1,55 +0,0 @@ -/* -*- c++ -*- ---------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifdef PAIR_CLASS - -PairStyle(gran/dmt/rolling,PairGranDMTRolling) - -#else - -#ifndef LMP_PAIR_GRAN_DMT_ROLLING_H -#define LMP_PAIR_GRAN_DMT_ROLLING_H - -#include "pair_gran_hooke_history.h" - -namespace LAMMPS_NS { - -class PairGranDMTRolling : public PairGranHookeHistory { -public: - PairGranDMTRolling(class LAMMPS *); - virtual ~PairGranDMTRolling(); - virtual void compute(int, int); - void settings(int, char **); //Eventually set this through coeff method so that user can specify a particular i-j set of coefficients - double single(int, int, int, int, double, double, double, double &); - double *E_one, *G_one, *pois, *muS_one, *cor, *alpha_one, *Ecoh_one, *kR_one, *muR_one, *etaR_one; //Public so as to be accessible to fix/wall/gran -private: - double **E, **G, **alpha, **muS, **Ecoh, **kR, **muR, **etaR, **gamman; - int normaldamp, rollingdamp; - - -}; - -} - -#endif -#endif - -/* ERROR/WARNING messages: - -E: Illegal ... command - -Self-explanatory. Check the input script syntax and compare to the -documentation for the command. You can use -echo screen as a -command-line option when running LAMMPS to see the offending line. - - */ diff --git a/src/GRANULAR/pair_gran_dmt_rolling2.cpp b/src/GRANULAR/pair_gran_dmt_rolling2.cpp deleted file mode 100644 index 5c1211cbc5..0000000000 --- a/src/GRANULAR/pair_gran_dmt_rolling2.cpp +++ /dev/null @@ -1,719 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. - ------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - Contributing authors: Leo Silbert (SNL), Gary Grest (SNL) - ------------------------------------------------------------------------- */ - -#include -#include -#include -#include -#include "pair_gran_dmt_rolling.h" -#include "atom.h" -#include "update.h" -#include "force.h" -#include "fix.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "comm.h" -#include "memory.h" -#include "error.h" -#include "math_const.h" - -using namespace LAMMPS_NS; -using namespace MathConst; - -#define TWOTHIRDS 0.6666666666666666 -#define EPSILON 1e-10 - -enum {TSUJI, BRILLIANTOV}; -enum {INDEP, BRILLROLL}; - -/* ---------------------------------------------------------------------- */ - -PairGranDMTRolling::PairGranDMTRolling(LAMMPS *lmp) : - PairGranHookeHistory(lmp, 7), - E_one(0), G_one(0), pois(0), muS_one(0), cor(0), alpha_one(0), - Ecoh_one(0), kR_one(0), muR_one(0), etaR_one(0) -{ - int ntypes = atom->ntypes; - memory->create(E,ntypes+1,ntypes+1,"pair:E"); - memory->create(G,ntypes+1,ntypes+1,"pair:G"); - memory->create(alpha,ntypes+1,ntypes+1,"pair:alpha"); - memory->create(gamman,ntypes+1,ntypes+1,"pair:gamman"); - memory->create(muS,ntypes+1,ntypes+1,"pair:muS"); - memory->create(Ecoh,ntypes+1,ntypes+1,"pair:Ecoh"); - memory->create(kR,ntypes+1,ntypes+1,"pair:kR"); - memory->create(muR,ntypes+1,ntypes+1,"pair:muR"); - memory->create(etaR,ntypes+1,ntypes+1,"pair:etaR"); -} - -/* ---------------------------------------------------------------------- */ -PairGranDMTRolling::~PairGranDMTRolling() -{ - delete [] E_one; - delete [] G_one; - delete [] pois; - delete [] muS_one; - delete [] cor; - delete [] alpha_one; - delete [] Ecoh_one; - delete [] kR_one; - delete [] muR_one; - delete [] etaR_one; - //TODO: Make all this work with standard pair coeff type commands. - //Also these should not be in the destructor. - memory->destroy(E); - memory->destroy(G); - memory->destroy(alpha); - memory->destroy(gamman); - memory->destroy(muS); - memory->destroy(Ecoh); - memory->destroy(kR); - memory->destroy(muR); - memory->destroy(etaR); -} -/* ---------------------------------------------------------------------- */ - -void PairGranDMTRolling::compute(int eflag, int vflag) -{ - int i,j,ii,jj,inum,jnum; - int itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,nx,ny,nz; - double radi,radj,radsum,rsq,r,rinv,rsqinv,R,a; - double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; - double wr1,wr2,wr3; - double vtr1,vtr2,vtr3,vrel; - double kn, kt, k_Q, k_R, eta_N, eta_T, eta_Q, eta_R; - double Fhz, Fdamp, Fdmt, Fne, Fntot, Fscrit, Frcrit; - double overlap; - double mi,mj,meff,damp,ccel,tor1,tor2,tor3; - double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3,vrlmag,vrlmaginv; - double rollmag, rolldotn, scalefac; - double fr, fr1, fr2, fr3; - double signtwist, magtwist, magtortwist, Mtcrit; - double fs,fs1,fs2,fs3,roll1,roll2,roll3,torroll1,torroll2,torroll3; - double tortwist1, tortwist2, tortwist3; - double shrmag,rsht; - int *ilist,*jlist,*numneigh,**firstneigh; - int *touch,**firsttouch; - double *shear,*allshear,**firstshear; - - if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = vflag_fdotr = 0; - - int shearupdate = 1; - if (update->setupflag) shearupdate = 0; - - // update rigid body info for owned & ghost atoms if using FixRigid masses - // body[i] = which body atom I is in, -1 if none - // mass_body = mass of each rigid body - - if (fix_rigid && neighbor->ago == 0){ - int tmp; - int *body = (int *) fix_rigid->extract("body",tmp); - double *mass_body = (double *) fix_rigid->extract("masstotal",tmp); - if (atom->nmax > nmax) { - memory->destroy(mass_rigid); - nmax = atom->nmax; - memory->create(mass_rigid,nmax,"pair:mass_rigid"); - } - int nlocal = atom->nlocal; - for (i = 0; i < nlocal; i++) - if (body[i] >= 0) mass_rigid[i] = mass_body[body[i]]; - else mass_rigid[i] = 0.0; - comm->forward_comm_pair(this); - } - - double **x = atom->x; - double **v = atom->v; - double **f = atom->f; - double **omega = atom->omega; - double **torque = atom->torque; - double *radius = atom->radius; - double *rmass = atom->rmass; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - inum = list->inum; - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->firstneigh; - firsttouch = list->listhistory->firstneigh; - firstshear = list->listhistory->firstdouble; - - // loop over neighbors of my atoms - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - itype = type[i]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - radi = radius[i]; - touch = firsttouch[i]; - allshear = firstshear[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - jtype = type[j]; - j &= NEIGHMASK; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - radj = radius[j]; - radsum = radi + radj; - - if (rsq >= radsum*radsum){ - // unset non-touching neighbors - touch[jj] = 0; - shear = &allshear[nsheardim*jj]; - for (int k = 0; k < nsheardim; k++) - shear[k] = 0.0; - } else { - r = sqrt(rsq); - rinv = 1.0/r; - rsqinv = 1.0/rsq; - R = radi*radj/(radi+radj); - nx = delx*rinv; - ny = dely*rinv; - nz = delz*rinv; - - // relative translational velocity - - vr1 = v[i][0] - v[j][0]; - vr2 = v[i][1] - v[j][1]; - vr3 = v[i][2] - v[j][2]; - - // normal component - - vnnr = vr1*nx + vr2*ny + vr3*nz; //v_R . n - vn1 = nx*vnnr; - vn2 = ny*vnnr; - vn3 = nz*vnnr; - - // meff = effective mass of pair of particles - // if I or J part of rigid body, use body mass - // if I or J is frozen, meff is other particle - - mi = rmass[i]; - mj = rmass[j]; - if (fix_rigid) { - if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; - if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; - } - - meff = mi*mj / (mi+mj); - if (mask[i] & freeze_group_bit) meff = mj; - if (mask[j] & freeze_group_bit) meff = mi; - - //**************************************** - //Normal force = Hertzian contact + DMT + damping - //**************************************** - overlap = radsum - r; - a = sqrt(R*overlap); - kn = 4.0/3.0*E[itype][jtype]*a; - Fhz = kn*overlap; - - //Damping (based on Tsuji et al) - if (normaldamp == BRILLIANTOV) eta_N = a*meff*gamman[itype][jtype]; - else if (normaldamp == TSUJI) eta_N=alpha[itype][jtype]*sqrt(meff*kn); - - Fdamp = -eta_N*vnnr; //F_nd eq 23 and Zhao eq 19 - - //DMT - Fdmt = -4*MY_PI*Ecoh[itype][jtype]*R; - - Fne = Fhz + Fdmt; - Fntot = Fne + Fdamp; - - //**************************************** - //Tangential force, including shear history effects - //**************************************** - - // tangential component - vt1 = vr1 - vn1; - vt2 = vr2 - vn2; - vt3 = vr3 - vn3; - - // relative rotational velocity - //Luding Gran Matt 2008, v10,p235 suggests correcting radi and radj by subtracting - //delta/2, i.e. instead of radi, use distance to center of contact point? - wr1 = (radi*omega[i][0] + radj*omega[j][0]); - wr2 = (radi*omega[i][1] + radj*omega[j][1]); - wr3 = (radi*omega[i][2] + radj*omega[j][2]); - - // relative tangential velocities - vtr1 = vt1 - (nz*wr2-ny*wr3); - vtr2 = vt2 - (nx*wr3-nz*wr1); - vtr3 = vt3 - (ny*wr1-nx*wr2); - vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; - vrel = sqrt(vrel); - - // shear history effects - touch[jj] = 1; - shear = &allshear[nsheardim*jj]; - shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + - shear[2]*shear[2]); - - // Rotate and update shear displacements. - // See e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 - if (shearupdate) { - rsht = shear[0]*nx + shear[1]*ny + shear[2]*nz; - if (fabs(rsht) < EPSILON) rsht = 0; - if (rsht > 0){ - scalefac = shrmag/(shrmag - rsht); //if rhst == shrmag, contacting pair has rotated 90 deg. in one step, in which case you deserve a crash! - shear[0] -= rsht*nx; - shear[1] -= rsht*ny; - shear[2] -= rsht*nz; - //Also rescale to preserve magnitude - shear[0] *= scalefac; - shear[1] *= scalefac; - shear[2] *= scalefac; - } - //Update shear history - shear[0] += vtr1*dt; - shear[1] += vtr2*dt; - shear[2] += vtr3*dt; - } - - // tangential forces = shear + tangential velocity damping - // following Zhao and Marshall Phys Fluids v20, p043302 (2008) - kt=8.0*G[itype][jtype]*a; - - eta_T = eta_N; //Based on discussion in Marshall; eta_T can also be an independent parameter - fs1 = -kt*shear[0] - eta_T*vtr1; //eq 26 - fs2 = -kt*shear[1] - eta_T*vtr2; - fs3 = -kt*shear[2] - eta_T*vtr3; - - // rescale frictional displacements and forces if needed - Fscrit = muS[itype][jtype] * fabs(Fne); - // For JKR, use eq 43 of Marshall. For DMT, use Fne instead - shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + - shear[2]*shear[2]); - fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); - if (fs > Fscrit) { - if (shrmag != 0.0) { - //shear[0] = (Fcrit/fs) * (shear[0] + eta_T*vtr1/kt) - eta_T*vtr1/kt; - //shear[1] = (Fcrit/fs) * (shear[1] + eta_T*vtr1/kt) - eta_T*vtr1/kt; - //shear[2] = (Fcrit/fs) * (shear[2] + eta_T*vtr1/kt) - eta_T*vtr1/kt; - shear[0] = -1.0/kt*(Fscrit*fs1/fs + eta_T*vtr1); //Same as above, but simpler (check!) - shear[1] = -1.0/kt*(Fscrit*fs2/fs + eta_T*vtr2); - shear[2] = -1.0/kt*(Fscrit*fs3/fs + eta_T*vtr3); - fs1 *= Fscrit/fs; - fs2 *= Fscrit/fs; - fs3 *= Fscrit/fs; - } else fs1 = fs2 = fs3 = 0.0; - } - - //**************************************** - // Rolling force, including shear history effects - //**************************************** - - relrot1 = omega[i][0] - omega[j][0]; - relrot2 = omega[i][1] - omega[j][1]; - relrot3 = omega[i][2] - omega[j][2]; - - // rolling velocity, see eq. 31 of Wang et al, Particuology v 23, p 49 (2015) - // This is different from the Marshall papers, which use the Bagi/Kuhn formulation - // for rolling velocity (see Wang et al for why the latter is wrong) - vrl1 = R*(relrot2*nz - relrot3*ny); //- 0.5*((radj-radi)/radsum)*vtr1; - vrl2 = R*(relrot3*nx - relrot1*nz); //- 0.5*((radj-radi)/radsum)*vtr2; - vrl3 = R*(relrot1*ny - relrot2*nx); //- 0.5*((radj-radi)/radsum)*vtr3; - vrlmag = sqrt(vrl1*vrl1+vrl2*vrl2+vrl3*vrl3); - if (vrlmag != 0.0) vrlmaginv = 1.0/vrlmag; - else vrlmaginv = 0.0; - - // Rolling displacement - rollmag = sqrt(shear[3]*shear[3] + shear[4]*shear[4] + shear[5]*shear[5]); - rolldotn = shear[3]*nx + shear[4]*ny + shear[5]*nz; - - if (shearupdate) { - if (fabs(rolldotn) < EPSILON) rolldotn = 0; - if (rolldotn > 0){ //Rotate into tangential plane - scalefac = rollmag/(rollmag - rolldotn); - shear[3] -= rolldotn*nx; - shear[4] -= rolldotn*ny; - shear[5] -= rolldotn*nz; - //Also rescale to preserve magnitude - shear[3] *= scalefac; - shear[4] *= scalefac; - shear[5] *= scalefac; - } - shear[3] += vrl1*dt; - shear[4] += vrl2*dt; - shear[5] += vrl3*dt; - } - - k_R = kR[itype][jtype]; - if (rollingdamp == INDEP) eta_R = etaR[itype][jtype]; - else if (rollingdamp == BRILLROLL) eta_R = muR[itype][jtype]*fabs(Fne); - fr1 = -k_R*shear[3] - eta_R*vrl1; - fr2 = -k_R*shear[4] - eta_R*vrl2; - fr3 = -k_R*shear[5] - eta_R*vrl3; - - // rescale frictional displacements and forces if needed - Frcrit = muR[itype][jtype] * fabs(Fne); - - fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3); - if (fr > Frcrit) { - if (rollmag != 0.0) { - shear[3] = -1.0/k_R*(Frcrit*fr1/fr + eta_R*vrl1); - shear[4] = -1.0/k_R*(Frcrit*fr2/fr + eta_R*vrl2); - shear[5] = -1.0/k_R*(Frcrit*fr3/fr + eta_R*vrl3); - fr1 *= Frcrit/fr; - fr2 *= Frcrit/fr; - fr3 *= Frcrit/fr; - } else fr1 = fr2 = fr3 = 0.0; - } - - - //**************************************** - // Twisting torque, including shear history effects - //**************************************** - magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of Marshall) - shear[6] += magtwist*dt; - k_Q = 0.5*kt*a*a;; //eq 32 - eta_Q = 0.5*eta_T*a*a; - magtortwist = -k_Q*shear[6] - eta_Q*magtwist;//M_t torque (eq 30) - - signtwist = (magtwist > 0) - (magtwist < 0); - Mtcrit=TWOTHIRDS*a*Fscrit;//critical torque (eq 44) - if (fabs(magtortwist) > Mtcrit){ - shear[6] = 1.0/k_Q*(Mtcrit*signtwist - eta_Q*magtwist); - magtortwist = -Mtcrit * signtwist; //eq 34 - } - - // Apply forces & torques - - fx = nx*Fntot + fs1; - fy = ny*Fntot + fs2; - fz = nz*Fntot + fs3; - - f[i][0] += fx; - f[i][1] += fy; - f[i][2] += fz; - - tor1 = ny*fs3 - nz*fs2; - tor2 = nz*fs1 - nx*fs3; - tor3 = nx*fs2 - ny*fs1; - - torque[i][0] -= radi*tor1; - torque[i][1] -= radi*tor2; - torque[i][2] -= radi*tor3; - - tortwist1 = magtortwist * nx; - tortwist2 = magtortwist * ny; - tortwist3 = magtortwist * nz; - - torque[i][0] += tortwist1; - torque[i][1] += tortwist2; - torque[i][2] += tortwist3; - - torroll1 = R*(ny*fr3 - nz*fr2); //n cross fr - torroll2 = R*(nz*fr1 - nx*fr3); - torroll3 = R*(nx*fr2 - ny*fr1); - - torque[i][0] += torroll1; - torque[i][1] += torroll2; - torque[i][2] += torroll3; - - if (force->newton_pair || j < nlocal) { - f[j][0] -= fx; - f[j][1] -= fy; - f[j][2] -= fz; - - torque[j][0] -= radj*tor1; - torque[j][1] -= radj*tor2; - torque[j][2] -= radj*tor3; - - torque[j][0] -= tortwist1; - torque[j][1] -= tortwist2; - torque[j][2] -= tortwist3; - - torque[j][0] -= torroll1; - torque[j][1] -= torroll2; - torque[j][2] -= torroll3; - } - if (evflag) ev_tally_xyz(i,j,nlocal,0, - 0.0,0.0,fx,fy,fz,delx,dely,delz); - } - } - } -} - -/* ---------------------------------------------------------------------- - global settings - ------------------------------------------------------------------------- */ - -void PairGranDMTRolling::settings(int narg, char **arg) -{ - if (narg < 6) error->all(FLERR,"Illegal pair_style command"); - - int ntypes = atom->ntypes; - - if (narg < 8*ntypes) error->all(FLERR,"Illegal pair_style command"); - - E_one = new double[ntypes+1]; - G_one = new double[ntypes+1]; - pois = new double[ntypes+1]; - muS_one = new double[ntypes+1]; - cor = new double[ntypes+1]; - alpha_one = new double[ntypes+1]; - Ecoh_one = new double[ntypes+1]; - kR_one = new double[ntypes+1]; - muR_one = new double[ntypes+1]; - etaR_one = new double[ntypes+1]; - - for (int i=0; i < ntypes;i++){ - E_one[i+1] = force->numeric(FLERR, arg[i]); - G_one[i+1] = force->numeric(FLERR, arg[ntypes+i]); - muS_one[i+1] = force->numeric(FLERR, arg[2*ntypes+i]); - cor[i+1] = force->numeric(FLERR, arg[3*ntypes+i]); - Ecoh_one[i+1] = force->numeric(FLERR, arg[4*ntypes+i]); - kR_one[i+1] = force->numeric(FLERR, arg[5*ntypes+i]); - muR_one[i+1] = force->numeric(FLERR, arg[6*ntypes+i]); - etaR_one[i+1] = force->numeric(FLERR, arg[7*ntypes+i]); - } - - //Defaults - normaldamp = TSUJI; - rollingdamp = INDEP; - - int iarg = 8*ntypes; - while (iarg < narg){ - if (strcmp(arg[iarg],"normaldamp") == 0){ - if (iarg+2 > narg) error->all(FLERR, "Invalid pair/gran/dmt/rolling entry"); - if (strcmp(arg[iarg+1],"tsuji") == 0) normaldamp = TSUJI; - else if (strcmp(arg[iarg+1],"brilliantov") == 0) normaldamp = BRILLIANTOV; - else error->all(FLERR, "Invalid normal damping model for pair/gran/dmt/rolling"); - iarg += 2; - } - else if (strcmp(arg[iarg],"rollingdamp") == 0){ - if (iarg+2 > narg) error->all(FLERR, "Invalid pair/gran/dmt/rolling entry"); - if (strcmp(arg[iarg+1],"independent") == 0) rollingdamp = INDEP; - else if (strcmp(arg[iarg+1],"brilliantov") == 0) rollingdamp = BRILLROLL; - else error->all(FLERR, "Invalid rolling damping model for pair/gran/dmt/rolling"); - iarg += 2; - } - else{ - iarg +=1; - } - } - - //Derived from inputs - for (int i=1; i <= ntypes; i++){ - pois[i] = E_one[i]/(2.0*G_one[i]) - 1.0; - alpha_one[i] = 1.2728-4.2783*cor[i]+11.087*cor[i]*cor[i]-22.348*cor[i]*cor[i]*cor[i]+27.467*cor[i]*cor[i]*cor[i]*cor[i]-18.022*cor[i]*cor[i]*cor[i]*cor[i]*cor[i]+4.8218*cor[i]*cor[i]*cor[i]*cor[i]*cor[i]*cor[i]; - for (int j=i; j <= ntypes; j++){ - E[i][j] = E[j][i] = 1/((1-pois[i]*pois[i])/E_one[i]+(1-pois[j]*pois[j])/E_one[j]); - G[i][j] = G[j][i] = 1/((2-pois[i])/G_one[i]+(2-pois[j])/G_one[j]); - if (normaldamp == TSUJI){ - alpha[i][j] = alpha[j][i] = sqrt(alpha_one[i]*alpha_one[j]); - } - else if (normaldamp == BRILLIANTOV){ - gamman[i][j] = gamman[j][i] = sqrt(cor[i]*cor[j]); - } - muS[i][j] = muS[j][i] = sqrt(muS_one[i]*muS_one[j]); - Ecoh[i][j] = Ecoh[j][i] = sqrt(Ecoh_one[i]*Ecoh_one[j]); - kR[i][j] = kR[j][i] = sqrt(kR_one[i]*kR_one[j]); - etaR[i][j] = etaR[j][i] = sqrt(etaR_one[i]*etaR_one[j]); - muR[i][j] = muR[j][i] = sqrt(muR_one[i]*muR_one[j]); - } - } -} - -/* ---------------------------------------------------------------------- */ - -double PairGranDMTRolling::single(int i, int j, int itype, int jtype, - double rsq, - double factor_coul, double factor_lj, - double &fforce) -{ - double radi,radj,radsum; - double r,rinv,rsqinv,delx,dely,delz, nx, ny, nz, R; - double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3,wr1,wr2,wr3; - double overlap, a; - double mi,mj,meff,damp,kn,kt; - double Fhz,Fdamp,Fdmt,Fne,Fntot,Fscrit; - double eta_N,eta_T; - double vtr1,vtr2,vtr3,vrel; - double fs1,fs2,fs3,fs; - double shrmag; - - - double *radius = atom->radius; - radi = radius[i]; - radj = radius[j]; - radsum = radi + radj; - - if (rsq >= radsum*radsum) { - fforce = 0.0; - svector[0] = svector[1] = svector[2] = svector[3] = 0.0; - return 0.0; - } - - r = sqrt(rsq); - rinv = 1.0/r; - rsqinv = 1.0/rsq; - R = radi*radj/radsum; - - // relative translational velocity - - double **v = atom->v; - vr1 = v[i][0] - v[j][0]; - vr2 = v[i][1] - v[j][1]; - vr3 = v[i][2] - v[j][2]; - - // normal component - - double **x = atom->x; - delx = x[i][0] - x[j][0]; - dely = x[i][1] - x[j][1]; - delz = x[i][2] - x[j][2]; - - nx = delx*rinv; - ny = dely*rinv; - nz = delz*rinv; - - - vnnr = vr1*nx + vr2*ny + vr3*nz; - vn1 = nx*vnnr; - vn2 = ny*vnnr; - vn3 = nz*vnnr; - - // tangential component - - vt1 = vr1 - vn1; - vt2 = vr2 - vn2; - vt3 = vr3 - vn3; - - // relative rotational velocity - - double **omega = atom->omega; - wr1 = (radi*omega[i][0] + radj*omega[j][0]); - wr2 = (radi*omega[i][1] + radj*omega[j][1]); - wr3 = (radi*omega[i][2] + radj*omega[j][2]); - - // meff = effective mass of pair of particles - // if I or J part of rigid body, use body mass - // if I or J is frozen, meff is other particle - - double *rmass = atom->rmass; - int *type = atom->type; - int *mask = atom->mask; - - mi = rmass[i]; - mj = rmass[j]; - if (fix_rigid) { - // NOTE: ensure mass_rigid is current for owned+ghost atoms? - if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; - if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; - } - - meff = mi*mj / (mi+mj); - if (mask[i] & freeze_group_bit) meff = mj; - if (mask[j] & freeze_group_bit) meff = mi; - - - // normal force = Hertzian contact + normal velocity damping - overlap = radsum - r; - a = sqrt(R*overlap); - kn = 4.0/3.0*E[itype][jtype]*a; - Fhz = kn*overlap; - - //Damping (based on Tsuji et al) - - eta_N=alpha[itype][jtype]*sqrt(meff*kn); - Fdamp = -eta_N*vnnr; //F_nd eq 23 and Zhao eq 19 - - //DMT - Fdmt = -4*MY_PI*Ecoh[itype][jtype]*R; - - Fne = Fhz + Fdmt; - Fntot = Fne + Fdamp; - - // relative velocities - - vtr1 = vt1 - (nz*wr2-ny*wr3); - vtr2 = vt2 - (nx*wr3-nz*wr1); - vtr3 = vt3 - (ny*wr1-nx*wr2); - vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; - vrel = sqrt(vrel); - - // shear history effects - // neighprev = index of found neigh on previous call - // search entire jnum list of neighbors of I for neighbor J - // start from neighprev, since will typically be next neighbor - // reset neighprev to 0 as necessary - - int jnum = list->numneigh[i]; - int *jlist = list->firstneigh[i]; - double *allshear = list->listhistory->firstdouble[i]; - - for (int jj = 0; jj < jnum; jj++) { - neighprev++; - if (neighprev >= jnum) neighprev = 0; - if (jlist[neighprev] == j) break; - } - - double *shear = &allshear[nsheardim*neighprev]; - shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + - shear[2]*shear[2]); - - // tangential forces = shear + tangential velocity damping - kt=8.0*G[itype][jtype]*a; - - eta_T = eta_N; - fs1 = -kt*shear[0] - eta_T*vtr1; - fs2 = -kt*shear[1] - eta_T*vtr2; - fs3 = -kt*shear[2] - eta_T*vtr3; - - // rescale frictional displacements and forces if needed - - fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); - Fscrit= muS[itype][jtype] * fabs(Fne); - - if (fs > Fscrit) { - if (shrmag != 0.0) { - fs1 *= Fscrit/fs; - fs2 *= Fscrit/fs; - fs3 *= Fscrit/fs; - fs *= Fscrit/fs; - } else fs1 = fs2 = fs3 = fs = 0.0; - } - - // set all forces and return no energy - - fforce = Fntot; - - // set single_extra quantities - - svector[0] = fs1; - svector[1] = fs2; - svector[2] = fs3; - svector[3] = fs; - svector[4] = vn1; - svector[5] = vn2; - svector[6] = vn3; - svector[7] = vt1; - svector[8] = vt2; - svector[9] = vt3; - return 0.0; -} diff --git a/src/GRANULAR/pair_gran_hooke_history_multi.cpp b/src/GRANULAR/pair_gran_hooke_history_multi.cpp deleted file mode 100644 index 48e793bbb3..0000000000 --- a/src/GRANULAR/pair_gran_hooke_history_multi.cpp +++ /dev/null @@ -1,915 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - Contributing authors: Leo Silbert (SNL), Gary Grest (SNL) -------------------------------------------------------------------------- */ - -#include -#include -#include -#include -#include "pair_gran_hooke_history_multi.h" -#include "atom.h" -#include "atom_vec.h" -#include "domain.h" -#include "force.h" -#include "update.h" -#include "modify.h" -#include "fix.h" -#include "fix_neigh_history.h" -#include "comm.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "neigh_request.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; - -#define BIG 1.0e20 - -/* ---------------------------------------------------------------------- */ - -PairGranHookeHistoryMulti::PairGranHookeHistoryMulti(LAMMPS *lmp) : Pair(lmp) -{ - single_enable = 1; - no_virial_fdotr_compute = 1; - history = 1; - fix_history = NULL; - - single_extra = 10; - svector = new double[10]; - - neighprev = 0; - - nmax = 0; - mass_rigid = NULL; - - // set comm size needed by this Pair if used with fix rigid - - comm_forward = 1; -} - -/* ---------------------------------------------------------------------- */ - -PairGranHookeHistoryMulti::~PairGranHookeHistoryMulti() -{ - delete [] svector; - if (fix_history) modify->delete_fix("NEIGH_HISTORY"); - - if (allocated) { - memory->destroy(setflag); - memory->destroy(cutsq); - - memory->destroy(cut); - memory->destroy(kn); - memory->destroy(kt); - memory->destroy(gamman); - memory->destroy(gammat); - memory->destroy(xmu); - memory->destroy(dampflag); - - delete [] onerad_dynamic; - delete [] onerad_frozen; - delete [] maxrad_dynamic; - delete [] maxrad_frozen; - } - memory->destroy(mass_rigid); -} - -/* ---------------------------------------------------------------------- */ - -void PairGranHookeHistoryMulti::compute(int eflag, int vflag) -{ - int i,j,ii,jj,inum,jnum,itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz; - double radi,radj,radsum,rsq,r,rinv,rsqinv; - double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; - double wr1,wr2,wr3; - double vtr1,vtr2,vtr3,vrel; - double mi,mj,meff,damp,ccel,tor1,tor2,tor3; - double fn,fs,fs1,fs2,fs3; - double shrmag,rsht; - int *ilist,*jlist,*numneigh,**firstneigh; - int *touch,**firsttouch; - double *shear,*allshear,**firstshear; - - if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = vflag_fdotr = 0; - - int shearupdate = 1; - if (update->setupflag) shearupdate = 0; - - // update rigid body info for owned & ghost atoms if using FixRigid masses - // body[i] = which body atom I is in, -1 if none - // mass_body = mass of each rigid body - - if (fix_rigid && neighbor->ago == 0) { - int tmp; - int *body = (int *) fix_rigid->extract("body",tmp); - double *mass_body = (double *) fix_rigid->extract("masstotal",tmp); - if (atom->nmax > nmax) { - memory->destroy(mass_rigid); - nmax = atom->nmax; - memory->create(mass_rigid,nmax,"pair:mass_rigid"); - } - int nlocal = atom->nlocal; - for (i = 0; i < nlocal; i++) - if (body[i] >= 0) mass_rigid[i] = mass_body[body[i]]; - else mass_rigid[i] = 0.0; - comm->forward_comm_pair(this); - } - - double **x = atom->x; - double **v = atom->v; - double **f = atom->f; - int *type = atom->type; - double **omega = atom->omega; - double **torque = atom->torque; - double *radius = atom->radius; - double *rmass = atom->rmass; - int *mask = atom->mask; - int nlocal = atom->nlocal; - int newton_pair = force->newton_pair; - - inum = list->inum; - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->firstneigh; - firsttouch = fix_history->firstflag; - firstshear = fix_history->firstvalue; - - // loop over neighbors of my atoms - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - itype = type[i]; - radi = radius[i]; - touch = firsttouch[i]; - allshear = firstshear[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - j &= NEIGHMASK; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - jtype = type[j]; - radj = radius[j]; - radsum = radi + radj; - - if (rsq >= radsum*radsum) { - - // unset non-touching neighbors - - touch[jj] = 0; - shear = &allshear[3*jj]; - shear[0] = 0.0; - shear[1] = 0.0; - shear[2] = 0.0; - - } else { - r = sqrt(rsq); - rinv = 1.0/r; - rsqinv = 1.0/rsq; - - // relative translational velocity - - vr1 = v[i][0] - v[j][0]; - vr2 = v[i][1] - v[j][1]; - vr3 = v[i][2] - v[j][2]; - - // normal component - - vnnr = vr1*delx + vr2*dely + vr3*delz; - vn1 = delx*vnnr * rsqinv; - vn2 = dely*vnnr * rsqinv; - vn3 = delz*vnnr * rsqinv; - - // tangential component - - vt1 = vr1 - vn1; - vt2 = vr2 - vn2; - vt3 = vr3 - vn3; - - // relative rotational velocity - - wr1 = (radi*omega[i][0] + radj*omega[j][0]) * rinv; - wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv; - wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv; - - // meff = effective mass of pair of particles - // if I or J part of rigid body, use body mass - // if I or J is frozen, meff is other particle - - mi = rmass[i]; - mj = rmass[j]; - if (fix_rigid) { - if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; - if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; - } - - meff = mi*mj / (mi+mj); - if (mask[i] & freeze_group_bit) meff = mj; - if (mask[j] & freeze_group_bit) meff = mi; - - // normal forces = Hookian contact + normal velocity damping - - damp = meff*gamman[itype][jtype]*vnnr*rsqinv; - ccel = kn[itype][jtype]*(radsum-r)*rinv - damp; - - // relative velocities - - vtr1 = vt1 - (delz*wr2-dely*wr3); - vtr2 = vt2 - (delx*wr3-delz*wr1); - vtr3 = vt3 - (dely*wr1-delx*wr2); - vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; - vrel = sqrt(vrel); - - // shear history effects - - touch[jj] = 1; - shear = &allshear[3*jj]; - - if (shearupdate) { - shear[0] += vtr1*dt; - shear[1] += vtr2*dt; - shear[2] += vtr3*dt; - } - shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + - shear[2]*shear[2]); - - // rotate shear displacements - - rsht = shear[0]*delx + shear[1]*dely + shear[2]*delz; - rsht *= rsqinv; - if (shearupdate) { - shear[0] -= rsht*delx; - shear[1] -= rsht*dely; - shear[2] -= rsht*delz; - } - - // tangential forces = shear + tangential velocity damping - - fs1 = - (kt[itype][jtype]*shear[0] + meff*gammat[itype][jtype]*vtr1); - fs2 = - (kt[itype][jtype]*shear[1] + meff*gammat[itype][jtype]*vtr2); - fs3 = - (kt[itype][jtype]*shear[2] + meff*gammat[itype][jtype]*vtr3); - - // rescale frictional displacements and forces if needed - - fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); - fn = xmu[itype][jtype] * fabs(ccel*r); - - if (fs > fn) { - if (shrmag != 0.0) { - shear[0] = (fn/fs) * (shear[0] + - meff*gammat[itype][jtype]*vtr1/kt[itype][jtype]) - - meff*gammat[itype][jtype]*vtr1/kt[itype][jtype]; - shear[1] = (fn/fs) * (shear[1] + - meff*gammat[itype][jtype]*vtr2/kt[itype][jtype]) - - meff*gammat[itype][jtype]*vtr2/kt[itype][jtype]; - shear[2] = (fn/fs) * (shear[2] + - meff*gammat[itype][jtype]*vtr3/kt[itype][jtype]) - - meff*gammat[itype][jtype]*vtr3/kt[itype][jtype]; - fs1 *= fn/fs; - fs2 *= fn/fs; - fs3 *= fn/fs; - } else fs1 = fs2 = fs3 = 0.0; - } - - // forces & torques - - fx = delx*ccel + fs1; - fy = dely*ccel + fs2; - fz = delz*ccel + fs3; - f[i][0] += fx; - f[i][1] += fy; - f[i][2] += fz; - - tor1 = rinv * (dely*fs3 - delz*fs2); - tor2 = rinv * (delz*fs1 - delx*fs3); - tor3 = rinv * (delx*fs2 - dely*fs1); - torque[i][0] -= radi*tor1; - torque[i][1] -= radi*tor2; - torque[i][2] -= radi*tor3; - - if (newton_pair || j < nlocal) { - f[j][0] -= fx; - f[j][1] -= fy; - f[j][2] -= fz; - torque[j][0] -= radj*tor1; - torque[j][1] -= radj*tor2; - torque[j][2] -= radj*tor3; - } - - if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair, - 0.0,0.0,fx,fy,fz,delx,dely,delz); - } - } - } - - if (vflag_fdotr) virial_fdotr_compute(); -} - -/* ---------------------------------------------------------------------- - allocate all arrays -------------------------------------------------------------------------- */ - -void PairGranHookeHistoryMulti::allocate() -{ - allocated = 1; - int n = atom->ntypes; - - memory->create(setflag,n+1,n+1,"pair:setflag"); - for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) - setflag[i][j] = 0; - - memory->create(cutsq,n+1,n+1,"pair:cutsq"); - memory->create(cut,n+1,n+1,"pair:cut"); - memory->create(kn,n+1,n+1,"pair:kn"); - memory->create(kt,n+1,n+1,"pair:kt"); - memory->create(gamman,n+1,n+1,"pair:gamman"); - memory->create(gammat,n+1,n+1,"pair:gammat"); - memory->create(xmu,n+1,n+1,"pair:xmu"); - memory->create(dampflag,n+1,n+1,"pair:dampflag"); - - onerad_dynamic = new double[n+1]; - onerad_frozen = new double[n+1]; - maxrad_dynamic = new double[n+1]; - maxrad_frozen = new double[n+1]; -} - -/* ---------------------------------------------------------------------- - global settings -------------------------------------------------------------------------- */ - -void PairGranHookeHistoryMulti::settings(int narg, char **arg) -{ - if (narg != 1) error->all(FLERR,"Illegal pair_style command"); - - if (strcmp(arg[0],"NULL") == 0 ) cut_global = -1.0; - else cut_global = force->numeric(FLERR,arg[0]); - - // reset cutoffs that have been explicitly set - if (allocated) { - int i,j; - for (i = 1; i <= atom->ntypes; i++) - for (j = i; j <= atom->ntypes; j++) - if (setflag[i][j]) cut[i][j] = cut_global; - } -} - -/* ---------------------------------------------------------------------- - set coeffs for one or more type pairs -------------------------------------------------------------------------- */ - -void PairGranHookeHistoryMulti::coeff(int narg, char **arg) -{ - if (narg < 8 || narg > 9) - error->all(FLERR,"Incorrect args for pair coefficients"); - - if (!allocated) allocate(); - - int ilo,ihi,jlo,jhi; - force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi); - force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi); - - double kn_one = force->numeric(FLERR,arg[2]); - double kt_one; - if (strcmp(arg[3],"NULL") == 0) kt_one = kn_one * 2.0/7.0; - else kt_one = force->numeric(FLERR,arg[3]); - - double gamman_one = force->numeric(FLERR,arg[4]); - double gammat_one; - if (strcmp(arg[5],"NULL") == 0) gammat_one = 0.5 * gamman_one; - else gammat_one = force->numeric(FLERR,arg[5]); - - double xmu_one = force->numeric(FLERR,arg[6]); - int dampflag_one = force->inumeric(FLERR,arg[7]); - if (dampflag_one == 0) gammat_one = 0.0; - - if (kn_one < 0.0 || kt_one < 0.0 || gamman_one < 0.0 || gammat_one < 0.0 || - xmu_one < 0.0 || xmu_one > 10000.0 || dampflag_one < 0 || dampflag_one > 1) - error->all(FLERR,"Illegal pair_style command"); - - // convert Kn and Kt from pressure units to force/distance^2 - kn_one /= force->nktv2p; - kt_one /= force->nktv2p; - - double cut_one = cut_global; - if (narg==9) { - if (strcmp(arg[8],"NULL") == 0) cut_one = -1.0; - else cut_one = force->numeric(FLERR,arg[8]); - } - - int count = 0; - for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { - kn[i][j] = kn_one; - kt[i][j] = kt_one; - gamman[i][j] = gamman_one; - gammat[i][j] = gammat_one; - xmu[i][j] = xmu_one; - dampflag[i][j] = dampflag_one; - cut[i][j] = cut_one; - setflag[i][j] = 1; - count++; - } - } - - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); -} - -/* ---------------------------------------------------------------------- - init specific to this pair style -------------------------------------------------------------------------- */ - -void PairGranHookeHistoryMulti::init_style() -{ - int i; - - // error and warning checks - - if (!atom->radius_flag || !atom->rmass_flag) - error->all(FLERR,"Pair granular requires atom attributes radius, rmass"); - if (comm->ghost_velocity == 0) - error->all(FLERR,"Pair granular requires ghost atoms store velocity"); - - // need a granular neigh list - - int irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->size = 1; - if (history) neighbor->requests[irequest]->history = 1; - - dt = update->dt; - - // if shear history is stored: - // if first init, create Fix needed for storing shear history - - if (history && fix_history == NULL) { - char dnumstr[16]; - sprintf(dnumstr,"%d",3); - char **fixarg = new char*[4]; - fixarg[0] = (char *) "NEIGH_HISTORY"; - fixarg[1] = (char *) "all"; - fixarg[2] = (char *) "NEIGH_HISTORY"; - fixarg[3] = dnumstr; - modify->add_fix(4,fixarg,1); - delete [] fixarg; - fix_history = (FixNeighHistory *) modify->fix[modify->nfix-1]; - fix_history->pair = this; - } - - // check for FixFreeze and set freeze_group_bit - - for (i = 0; i < modify->nfix; i++) - if (strcmp(modify->fix[i]->style,"freeze") == 0) break; - if (i < modify->nfix) freeze_group_bit = modify->fix[i]->groupbit; - else freeze_group_bit = 0; - - // check for FixRigid so can extract rigid body masses - - fix_rigid = NULL; - for (i = 0; i < modify->nfix; i++) - if (modify->fix[i]->rigid_flag) break; - if (i < modify->nfix) fix_rigid = modify->fix[i]; - - // check for FixPour and FixDeposit so can extract particle radii - - int ipour; - for (ipour = 0; ipour < modify->nfix; ipour++) - if (strcmp(modify->fix[ipour]->style,"pour") == 0) break; - if (ipour == modify->nfix) ipour = -1; - - int idep; - for (idep = 0; idep < modify->nfix; idep++) - if (strcmp(modify->fix[idep]->style,"deposit") == 0) break; - if (idep == modify->nfix) idep = -1; - - // set maxrad_dynamic and maxrad_frozen for each type - // include future FixPour and FixDeposit particles as dynamic - - int itype; - for (i = 1; i <= atom->ntypes; i++) { - onerad_dynamic[i] = onerad_frozen[i] = 0.0; - if (ipour >= 0) { - itype = i; - onerad_dynamic[i] = - *((double *) modify->fix[ipour]->extract("radius",itype)); - } - if (idep >= 0) { - itype = i; - onerad_dynamic[i] = - *((double *) modify->fix[idep]->extract("radius",itype)); - } - } - - double *radius = atom->radius; - int *mask = atom->mask; - int *type = atom->type; - int nlocal = atom->nlocal; - - for (i = 0; i < nlocal; i++) - if (mask[i] & freeze_group_bit) - onerad_frozen[type[i]] = MAX(onerad_frozen[type[i]],radius[i]); - else - onerad_dynamic[type[i]] = MAX(onerad_dynamic[type[i]],radius[i]); - - MPI_Allreduce(&onerad_dynamic[1],&maxrad_dynamic[1],atom->ntypes, - MPI_DOUBLE,MPI_MAX,world); - MPI_Allreduce(&onerad_frozen[1],&maxrad_frozen[1],atom->ntypes, - MPI_DOUBLE,MPI_MAX,world); - - // set fix which stores history info - - if (history) { - int ifix = modify->find_fix("NEIGH_HISTORY"); - if (ifix < 0) error->all(FLERR,"Could not find pair fix neigh history ID"); - fix_history = (FixNeighHistory *) modify->fix[ifix]; - } -} - -/* ---------------------------------------------------------------------- - init for one type pair i,j and corresponding j,i -------------------------------------------------------------------------- */ - -double PairGranHookeHistoryMulti::init_one(int i, int j) -{ - if (setflag[i][j] == 0) { - kn[i][j] = mix_stiffness(kn[i][i],kn[j][j]); - kt[i][j] = mix_stiffness(kt[i][i],kt[j][j]); - gamman[i][j] = mix_damping(gamman[i][i],gamman[j][j]); - gammat[i][j] = mix_damping(gammat[i][i],gammat[j][j]); - xmu[i][j] = mix_friction(xmu[i][i],xmu[j][j]); - - dampflag[i][j] = 0; - if (dampflag[i][i] || dampflag[j][j]) dampflag[i][j] = 1; - - } - - kn[j][i] = kn[i][j]; - kt[j][i] = kt[i][j]; - gamman[j][i] = gamman[i][j]; - gammat[j][i] = gammat[i][j]; - xmu[j][i] = xmu[i][j]; - dampflag[j][i] = dampflag[i][j]; - - double cutoff = cut[i][j]; - - // It is likely that cut[i][j] at this point is still 0.0. This can happen when - // there is a future fix_pour after the current run. A cut[i][j] = 0.0 creates - // problems because neighbor.cpp uses min(cut[i][j]) to decide on the bin size - // To avoid this issue,for cases involving cut[i][j] = 0.0 (possible only - // if there is no current information about radius/cutoff of type i and j). - // we assign cutoff = min(cut[i][j]) for i,j such that cut[i][j] > 0.0. - - if (cut[i][j] < 0.0) { - if (((maxrad_dynamic[i] > 0.0) && (maxrad_dynamic[j] > 0.0)) || ((maxrad_dynamic[i] > 0.0) && (maxrad_frozen[j] > 0.0)) || - ((maxrad_frozen[i] > 0.0) && (maxrad_dynamic[j] > 0.0))) { // radius info about both i and j exist - cutoff = maxrad_dynamic[i]+maxrad_dynamic[j]; - cutoff = MAX(cutoff,maxrad_frozen[i]+maxrad_dynamic[j]); - cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]); - } - else { // radius info about either i or j does not exist (i.e. not present and not about to get poured; set to largest value to not interfere with neighbor list) - double cutmax = 0.0; - for (int k = 1; k <= atom->ntypes; k++) { - cutmax = MAX(cutmax,2.0*maxrad_dynamic[k]); - cutmax = MAX(cutmax,2.0*maxrad_frozen[k]); - } - cutoff = cutmax; - } - } - return cutoff; -} - -/* ---------------------------------------------------------------------- - proc 0 writes to restart file -------------------------------------------------------------------------- */ - -void PairGranHookeHistoryMulti::write_restart(FILE *fp) -{ - write_restart_settings(fp); - - int i,j; - for (i = 1; i <= atom->ntypes; i++) { - for (j = i; j <= atom->ntypes; j++) { - fwrite(&setflag[i][j],sizeof(int),1,fp); - if (setflag[i][j]) { - fwrite(&kn[i][j],sizeof(double),1,fp); - fwrite(&kt[i][j],sizeof(double),1,fp); - fwrite(&gamman[i][j],sizeof(double),1,fp); - fwrite(&gammat[i][j],sizeof(double),1,fp); - fwrite(&xmu[i][j],sizeof(double),1,fp); - fwrite(&dampflag[i][j],sizeof(int),1,fp); - fwrite(&cut[i][j],sizeof(double),1,fp); - } - } - } -} - -/* ---------------------------------------------------------------------- - proc 0 reads from restart file, bcasts -------------------------------------------------------------------------- */ - -void PairGranHookeHistoryMulti::read_restart(FILE *fp) -{ - read_restart_settings(fp); - allocate(); - - int i,j; - int me = comm->me; - for (i = 1; i <= atom->ntypes; i++) { - for (j = i; j <= atom->ntypes; j++) { - if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); - MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); - if (setflag[i][j]) { - if (me == 0) { - fread(&kn[i][j],sizeof(double),1,fp); - fread(&kt[i][j],sizeof(double),1,fp); - fread(&gamman[i][j],sizeof(double),1,fp); - fread(&gammat[i][j],sizeof(double),1,fp); - fread(&xmu[i][j],sizeof(double),1,fp); - fread(&dampflag[i][j],sizeof(int),1,fp); - fread(&cut[i][j],sizeof(double),1,fp); - } - MPI_Bcast(&kn[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&kt[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&gamman[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&gammat[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&xmu[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&dampflag[i][j],1,MPI_INT,0,world); - } - } - } -} - -/* ---------------------------------------------------------------------- - proc 0 writes to restart file -------------------------------------------------------------------------- */ - -void PairGranHookeHistoryMulti::write_restart_settings(FILE *fp) -{ - fwrite(&cut_global,sizeof(double),1,fp); -} - -/* ---------------------------------------------------------------------- - proc 0 reads from restart file, bcasts -------------------------------------------------------------------------- */ - -void PairGranHookeHistoryMulti::read_restart_settings(FILE *fp) -{ - if (comm->me == 0) { - fread(&cut_global,sizeof(double),1,fp); - } - MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); -} - -/* ---------------------------------------------------------------------- */ - -void PairGranHookeHistoryMulti::reset_dt() -{ - dt = update->dt; -} - -/* ---------------------------------------------------------------------- */ - -double PairGranHookeHistoryMulti::single(int i, int j, int itype, int jtype, - double rsq, - double factor_coul, double factor_lj, - double &fforce) -{ - double radi,radj,radsum; - double r,rinv,rsqinv,delx,dely,delz; - double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3,wr1,wr2,wr3; - double mi,mj,meff,damp,ccel; - double vtr1,vtr2,vtr3,vrel,shrmag,rsht; - double fs1,fs2,fs3,fs,fn; - - double *radius = atom->radius; - radi = radius[i]; - radj = radius[j]; - radsum = radi + radj; - - if (rsq >= radsum*radsum) { - fforce = 0.0; - for (int m = 0; m < single_extra; m++) svector[m] = 0.0; - return 0.0; - } - - r = sqrt(rsq); - rinv = 1.0/r; - rsqinv = 1.0/rsq; - - // relative translational velocity - - double **v = atom->v; - vr1 = v[i][0] - v[j][0]; - vr2 = v[i][1] - v[j][1]; - vr3 = v[i][2] - v[j][2]; - - // normal component - - double **x = atom->x; - delx = x[i][0] - x[j][0]; - dely = x[i][1] - x[j][1]; - delz = x[i][2] - x[j][2]; - - vnnr = vr1*delx + vr2*dely + vr3*delz; - vn1 = delx*vnnr * rsqinv; - vn2 = dely*vnnr * rsqinv; - vn3 = delz*vnnr * rsqinv; - - // tangential component - - vt1 = vr1 - vn1; - vt2 = vr2 - vn2; - vt3 = vr3 - vn3; - - // relative rotational velocity - - double **omega = atom->omega; - wr1 = (radi*omega[i][0] + radj*omega[j][0]) * rinv; - wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv; - wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv; - - // meff = effective mass of pair of particles - // if I or J part of rigid body, use body mass - // if I or J is frozen, meff is other particle - - double *rmass = atom->rmass; - int *mask = atom->mask; - - mi = rmass[i]; - mj = rmass[j]; - if (fix_rigid) { - // NOTE: insure mass_rigid is current for owned+ghost atoms? - if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; - if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; - } - - meff = mi*mj / (mi+mj); - if (mask[i] & freeze_group_bit) meff = mj; - if (mask[j] & freeze_group_bit) meff = mi; - - // normal forces = Hookian contact + normal velocity damping - - damp = meff*gamman[itype][jtype]*vnnr*rsqinv; - ccel = kn[itype][jtype]*(radsum-r)*rinv - damp; - - // relative velocities - - vtr1 = vt1 - (delz*wr2-dely*wr3); - vtr2 = vt2 - (delx*wr3-delz*wr1); - vtr3 = vt3 - (dely*wr1-delx*wr2); - vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; - vrel = sqrt(vrel); - - // shear history effects - // neighprev = index of found neigh on previous call - // search entire jnum list of neighbors of I for neighbor J - // start from neighprev, since will typically be next neighbor - // reset neighprev to 0 as necessary - - int jnum = list->numneigh[i]; - int *jlist = list->firstneigh[i]; - double *allshear = fix_history->firstvalue[i]; - - for (int jj = 0; jj < jnum; jj++) { - neighprev++; - if (neighprev >= jnum) neighprev = 0; - if (jlist[neighprev] == j) break; - } - - double *shear = &allshear[3*neighprev]; - shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + - shear[2]*shear[2]); - - // rotate shear displacements - - rsht = shear[0]*delx + shear[1]*dely + shear[2]*delz; - rsht *= rsqinv; - - // tangential forces = shear + tangential velocity damping - - fs1 = - (kt[itype][jtype]*shear[0] + meff*gammat[itype][jtype]*vtr1); - fs2 = - (kt[itype][jtype]*shear[1] + meff*gammat[itype][jtype]*vtr2); - fs3 = - (kt[itype][jtype]*shear[2] + meff*gammat[itype][jtype]*vtr3); - - // rescale frictional displacements and forces if needed - - fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); - fn = xmu[itype][jtype] * fabs(ccel*r); - - if (fs > fn) { - if (shrmag != 0.0) { - fs1 *= fn/fs; - fs2 *= fn/fs; - fs3 *= fn/fs; - fs *= fn/fs; - } else fs1 = fs2 = fs3 = fs = 0.0; - } - - // set force and return no energy - - fforce = ccel; - - // set single_extra quantities - - svector[0] = fs1; - svector[1] = fs2; - svector[2] = fs3; - svector[3] = fs; - svector[4] = vn1; - svector[5] = vn2; - svector[6] = vn3; - svector[7] = vt1; - svector[8] = vt2; - svector[9] = vt3; - - return 0.0; -} - -/* ---------------------------------------------------------------------- */ - -int PairGranHookeHistoryMulti::pack_forward_comm(int n, int *list, double *buf, - int pbc_flag, int *pbc) -{ - int i,j,m; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = mass_rigid[j]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -void PairGranHookeHistoryMulti::unpack_forward_comm(int n, int first, double *buf) -{ - int i,m,last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) - mass_rigid[i] = buf[m++]; -} - -/* ---------------------------------------------------------------------- - memory usage of local atom-based arrays -------------------------------------------------------------------------- */ - -double PairGranHookeHistoryMulti::memory_usage() -{ - double bytes = nmax * sizeof(double); - return bytes; -} - -/* ---------------------------------------------------------------------- - mixing of stiffness -------------------------------------------------------------------------- */ - -double PairGranHookeHistoryMulti::mix_stiffness(double kii, double kjj) -{ - return kii*kjj/(kii + kjj); -} - -/* ---------------------------------------------------------------------- - mixing of damping -------------------------------------------------------------------------- */ - -double PairGranHookeHistoryMulti::mix_damping(double gammaii, double gammajj) -{ - return sqrt(gammaii*gammajj); -} - -/* ---------------------------------------------------------------------- - mixing of friction -------------------------------------------------------------------------- */ - -double PairGranHookeHistoryMulti::mix_friction(double xmuii, double xmujj) -{ - return MAX(xmuii,xmujj); -} - diff --git a/src/GRANULAR/pair_gran_hooke_history_multi.h b/src/GRANULAR/pair_gran_hooke_history_multi.h deleted file mode 100644 index f302ede96c..0000000000 --- a/src/GRANULAR/pair_gran_hooke_history_multi.h +++ /dev/null @@ -1,109 +0,0 @@ -/* -*- c++ -*- ---------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifdef PAIR_CLASS - -PairStyle(gran/hooke/history/multi,PairGranHookeHistoryMulti) - -#else - -#ifndef LMP_PAIR_GRAN_HOOKE_HISTORY_MULTI_H -#define LMP_PAIR_GRAN_HOOKE_HISTORY_MULTI_H - -#include "pair.h" - -namespace LAMMPS_NS { - -class PairGranHookeHistoryMulti : public Pair { - public: - PairGranHookeHistoryMulti(class LAMMPS *); - virtual ~PairGranHookeHistoryMulti(); - virtual void compute(int, int); - virtual void settings(int, char **); - virtual void coeff(int, char **); // Made Virtual by IS Oct 7 2017 - void init_style(); - double init_one(int, int); - void write_restart(FILE *); - void read_restart(FILE *); - void write_restart_settings(FILE *); - void read_restart_settings(FILE *); - void reset_dt(); - virtual double single(int, int, int, int, double, double, double, double &); - int pack_forward_comm(int, int *, double *, int, int *); - void unpack_forward_comm(int, int, double *); - double memory_usage(); - - protected: - double cut_global; - double **kn,**kt,**gamman,**gammat,**xmu,**cut; - int **dampflag; - double dt; - int freeze_group_bit; - int history; - - int neighprev; - double *onerad_dynamic,*onerad_frozen; - double *maxrad_dynamic,*maxrad_frozen; - - class FixNeighHistory *fix_history; - - // storage of rigid body masses for use in granular interactions - - class Fix *fix_rigid; // ptr to rigid body fix, NULL if none - double *mass_rigid; // rigid mass for owned+ghost atoms - int nmax; // allocated size of mass_rigid - - virtual void allocate(); // Made Virtual by IS Oct 7 2017 - -private: - double mix_stiffness(double kii, double kjj); - double mix_damping(double gammaii, double gammajj); - double mix_friction(double xmuii, double xmujj); -}; - -} - -#endif -#endif - -/* ERROR/WARNING messages: - -E: Illegal ... command - -Self-explanatory. Check the input script syntax and compare to the -documentation for the command. You can use -echo screen as a -command-line option when running LAMMPS to see the offending line. - -E: Incorrect args for pair coefficients - -Self-explanatory. Check the input script or data file. - -E: Pair granular requires atom attributes radius, rmass - -The atom style defined does not have these attributes. - -E: Pair granular requires ghost atoms store velocity - -Use the comm_modify vel yes command to enable this. - -E: Pair granular with shear history requires newton pair off - -This is a current restriction of the implementation of pair -granular styles with history. - -E: Could not find pair fix ID - -A fix is created internally by the pair style to store shear -history information. You cannot delete it. - -*/ diff --git a/src/GRANULAR/pair_gran_jkr_rolling.cpp b/src/GRANULAR/pair_gran_jkr_rolling.cpp deleted file mode 100644 index ce109cccbc..0000000000 --- a/src/GRANULAR/pair_gran_jkr_rolling.cpp +++ /dev/null @@ -1,763 +0,0 @@ -/* ---------------------------------------------------------------------- - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. - ------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - Contributing authors: Leo Silbert (SNL), Gary Grest (SNL) - ------------------------------------------------------------------------- */ - -#include -#include -#include -#include -#include "pair_gran_jkr_rolling.h" -#include "atom.h" -#include "update.h" -#include "force.h" -#include "fix.h" -#include "fix_neigh_history.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "comm.h" -#include "memory.h" -#include "error.h" -#include "math_const.h" - -using namespace LAMMPS_NS; -using namespace MathConst; - -#define ONETHIRD 0.33333333333333333 -#define TWOTHIRDS 0.66666666666666666 -#define POW6ONE 0.550321208149104 //6^(-1/3) -#define POW6TWO 0.30285343213869 //6^(-2/3) - -#define EPSILON 1e-10 - -enum {TSUJI, BRILLIANTOV}; -enum {INDEP, BRILLROLL}; - -/* ---------------------------------------------------------------------- */ - -PairGranJKRRolling::PairGranJKRRolling(LAMMPS *lmp) : - PairGranHookeHistory(lmp, 7) -{ - E_one = NULL; - G_one = NULL; - pois = NULL; - muS_one = NULL; - cor = NULL; - alpha_one = NULL; - Ecoh_one = NULL; - kR_one = NULL; - muR_one = NULL; - etaR_one = NULL; - int ntypes = atom->ntypes; - memory->create(E,ntypes+1,ntypes+1,"pair:E"); - memory->create(G,ntypes+1,ntypes+1,"pair:G"); - memory->create(alpha,ntypes+1,ntypes+1,"pair:alpha"); - memory->create(gamman,ntypes+1,ntypes+1,"pair:gamman"); - memory->create(muS,ntypes+1,ntypes+1,"pair:muS"); - memory->create(Ecoh,ntypes+1,ntypes+1,"pair:Ecoh"); - memory->create(kR,ntypes+1,ntypes+1,"pair:kR"); - memory->create(muR,ntypes+1,ntypes+1,"pair:muR"); - memory->create(etaR,ntypes+1,ntypes+1,"pair:etaR"); -} - -/* ---------------------------------------------------------------------- */ -PairGranJKRRolling::~PairGranJKRRolling() -{ - delete [] E_one; - delete [] G_one; - delete [] pois; - delete [] muS_one; - delete [] cor; - delete [] alpha_one; - delete [] Ecoh_one; - delete [] kR_one; - delete [] muR_one; - delete [] etaR_one; - //TODO: Make all this work with standard pair coeff type commands. - //Also these should not be in the destructor. - memory->destroy(E); - memory->destroy(G); - memory->destroy(alpha); - memory->destroy(gamman); - memory->destroy(muS); - memory->destroy(Ecoh); - memory->destroy(kR); - memory->destroy(muR); - memory->destroy(etaR); -} -/* ---------------------------------------------------------------------- */ - -void PairGranJKRRolling::compute(int eflag, int vflag) -{ - int i,j,ii,jj,inum,jnum; - int itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,nx,ny,nz; - double radi,radj,radsum,rsq,r,rinv,rsqinv,R,a; - double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; - double wr1,wr2,wr3; - double vtr1,vtr2,vtr3,vrel; - double kn, kt, k_Q, k_R, eta_N, eta_T, eta_Q, eta_R; - double Fne, Fdamp, Fntot, Fscrit, Frcrit, F_C, delta_C, delta_Cinv; - double overlap, olapsq, olapcubed, sqrtterm, tmp, a0; - double keyterm, keyterm2, keyterm3, aovera0, foverFc; - double mi,mj,meff,damp,ccel,tor1,tor2,tor3; - double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3,vrlmag,vrlmaginv; - double rollmag, rolldotn, scalefac; - double fr, fr1, fr2, fr3; - double signtwist, magtwist, magtortwist, Mtcrit; - double fs,fs1,fs2,fs3,roll1,roll2,roll3,torroll1,torroll2,torroll3; - double tortwist1, tortwist2, tortwist3; - double shrmag,rsht; - int *ilist,*jlist,*numneigh,**firstneigh; - int *touch,**firsttouch; - double *shear,*allshear,**firstshear; - - if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = vflag_fdotr = 0; - - int shearupdate = 1; - if (update->setupflag) shearupdate = 0; - - // update rigid body info for owned & ghost atoms if using FixRigid masses - // body[i] = which body atom I is in, -1 if none - // mass_body = mass of each rigid body - - if (fix_rigid && neighbor->ago == 0){ - int tmp; - int *body = (int *) fix_rigid->extract("body",tmp); - double *mass_body = (double *) fix_rigid->extract("masstotal",tmp); - if (atom->nmax > nmax) { - memory->destroy(mass_rigid); - nmax = atom->nmax; - memory->create(mass_rigid,nmax,"pair:mass_rigid"); - } - int nlocal = atom->nlocal; - for (i = 0; i < nlocal; i++) - if (body[i] >= 0) mass_rigid[i] = mass_body[body[i]]; - else mass_rigid[i] = 0.0; - comm->forward_comm_pair(this); - } - - double **x = atom->x; - double **v = atom->v; - double **f = atom->f; - double **omega = atom->omega; - double **torque = atom->torque; - double *radius = atom->radius; - double *rmass = atom->rmass; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - inum = list->inum; - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->firstneigh; - firsttouch = fix_history->firstflag; - firstshear = fix_history->firstvalue; - - // loop over neighbors of my atoms - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - itype = type[i]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - radi = radius[i]; - touch = firsttouch[i]; - allshear = firstshear[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - jtype = type[j]; - j &= NEIGHMASK; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - radj = radius[j]; - radsum = radi + radj; - R = radi*radj/(radi+radj); - a0 = pow(9.0*M_PI*Ecoh[itype][jtype]*R*R/E[itype][jtype],ONETHIRD); - delta_C = 0.5*a0*a0*POW6ONE/R; - - if ((rsq >= radsum*radsum && touch[jj] == 0) || - (rsq >= (radsum+delta_C)*(radsum+delta_C))){ - // unset non-touching neighbors - touch[jj] = 0; - shear = &allshear[size_history*jj]; - for (int k = 0; k < size_history; k++) - shear[k] = 0.0; - } else { - F_C = 3.0*R*M_PI*Ecoh[itype][jtype]; - r = sqrt(rsq); - rinv = 1.0/r; - rsqinv = 1.0/rsq; - - nx = delx*rinv; - ny = dely*rinv; - nz = delz*rinv; - - // relative translational velocity - - vr1 = v[i][0] - v[j][0]; - vr2 = v[i][1] - v[j][1]; - vr3 = v[i][2] - v[j][2]; - - // normal component - - vnnr = vr1*nx + vr2*ny + vr3*nz; //v_R . n - vn1 = nx*vnnr; - vn2 = ny*vnnr; - vn3 = nz*vnnr; - - // meff = effective mass of pair of particles - // if I or J part of rigid body, use body mass - // if I or J is frozen, meff is other particle - - mi = rmass[i]; - mj = rmass[j]; - if (fix_rigid) { - if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; - if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; - } - - meff = mi*mj / (mi+mj); - if (mask[i] & freeze_group_bit) meff = mj; - if (mask[j] & freeze_group_bit) meff = mi; - - //**************************************** - //Normal force = JKR-adjusted Hertzian contact + damping - //**************************************** - if (Ecoh[itype][jtype] != 0.0) delta_Cinv = 1.0/delta_C; - else delta_Cinv = 1.0; - overlap = (radsum - r)*delta_Cinv; - olapsq = overlap*overlap; - olapcubed = olapsq*overlap; - sqrtterm = sqrt(1.0 + olapcubed); - tmp = 2.0 + olapcubed + 2.0*sqrtterm; - keyterm = pow(tmp,ONETHIRD); - keyterm2 = olapsq/keyterm; - keyterm3 = sqrt(overlap + keyterm2 + keyterm); - aovera0 = POW6TWO * (keyterm3 + - sqrt(2.0*overlap - keyterm2 - keyterm + 4.0/keyterm3));// eq 41 - a = aovera0*a0; - foverFc = 4.0*((aovera0*aovera0*aovera0) - pow(aovera0,1.5));//F_ne/F_C (eq 40) - - Fne = F_C*foverFc; - - //Damping - kn = 4.0/3.0*E[itype][jtype]*a; - if (normaldamp == BRILLIANTOV) eta_N = a*meff*gamman[itype][jtype]; - else if (normaldamp == TSUJI) eta_N=alpha[itype][jtype]*sqrt(meff*kn); - - Fdamp = -eta_N*vnnr; //F_nd eq 23 and Zhao eq 19 - - Fntot = Fne + Fdamp; - - //**************************************** - //Tangential force, including shear history effects - //**************************************** - - // tangential component - vt1 = vr1 - vn1; - vt2 = vr2 - vn2; - vt3 = vr3 - vn3; - - // relative rotational velocity - //Luding Gran Matt 2008, v10,p235 suggests correcting radi and radj by subtracting - //delta/2, i.e. instead of radi, use distance to center of contact point? - wr1 = (radi*omega[i][0] + radj*omega[j][0]); - wr2 = (radi*omega[i][1] + radj*omega[j][1]); - wr3 = (radi*omega[i][2] + radj*omega[j][2]); - - // relative tangential velocities - vtr1 = vt1 - (nz*wr2-ny*wr3); - vtr2 = vt2 - (nx*wr3-nz*wr1); - vtr3 = vt3 - (ny*wr1-nx*wr2); - vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; - vrel = sqrt(vrel); - - // shear history effects - touch[jj] = 1; - shear = &allshear[size_history*jj]; - shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + - shear[2]*shear[2]); - - // Rotate and update shear displacements. - // See e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 - if (shearupdate) { - rsht = shear[0]*nx + shear[1]*ny + shear[2]*nz; - if (fabs(rsht) < EPSILON) rsht = 0; - if (rsht > 0){ - scalefac = shrmag/(shrmag - rsht); //if rhst == shrmag, contacting pair has rotated 90 deg. in one step, in which case you deserve a crash! - shear[0] -= rsht*nx; - shear[1] -= rsht*ny; - shear[2] -= rsht*nz; - //Also rescale to preserve magnitude - shear[0] *= scalefac; - shear[1] *= scalefac; - shear[2] *= scalefac; - } - //Update shear history - shear[0] += vtr1*dt; - shear[1] += vtr2*dt; - shear[2] += vtr3*dt; - } - - // tangential forces = shear + tangential velocity damping - // following Zhao and Marshall Phys Fluids v20, p043302 (2008) - kt=8.0*G[itype][jtype]*a; - - eta_T = eta_N; //Based on discussion in Marshall; eta_T can also be an independent parameter - fs1 = -kt*shear[0] - eta_T*vtr1; //eq 26 - fs2 = -kt*shear[1] - eta_T*vtr2; - fs3 = -kt*shear[2] - eta_T*vtr3; - - // rescale frictional displacements and forces if needed - Fscrit = muS[itype][jtype] * fabs(Fne + 2*F_C); - // For JKR, use eq 43 of Marshall. For DMT, use Fne instead - - fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); - if (fs > Fscrit) { - if (shrmag != 0.0) { - //shear[0] = (Fcrit/fs) * (shear[0] + eta_T*vtr1/kt) - eta_T*vtr1/kt; - //shear[1] = (Fcrit/fs) * (shear[1] + eta_T*vtr1/kt) - eta_T*vtr1/kt; - //shear[2] = (Fcrit/fs) * (shear[2] + eta_T*vtr1/kt) - eta_T*vtr1/kt; - shear[0] = -1.0/kt*(Fscrit*fs1/fs + eta_T*vtr1); //Same as above, but simpler (check!) - shear[1] = -1.0/kt*(Fscrit*fs2/fs + eta_T*vtr2); - shear[2] = -1.0/kt*(Fscrit*fs3/fs + eta_T*vtr3); - fs1 *= Fscrit/fs; - fs2 *= Fscrit/fs; - fs3 *= Fscrit/fs; - } else fs1 = fs2 = fs3 = 0.0; - } - - //**************************************** - // Rolling force, including shear history effects - //**************************************** - - relrot1 = omega[i][0] - omega[j][0]; - relrot2 = omega[i][1] - omega[j][1]; - relrot3 = omega[i][2] - omega[j][2]; - - // rolling velocity, see eq. 31 of Wang et al, Particuology v 23, p 49 (2015) - // This is different from the Marshall papers, which use the Bagi/Kuhn formulation - // for rolling velocity (see Wang et al for why the latter is wrong) - vrl1 = R*(relrot2*nz - relrot3*ny); //- 0.5*((radj-radi)/radsum)*vtr1; - vrl2 = R*(relrot3*nx - relrot1*nz); //- 0.5*((radj-radi)/radsum)*vtr2; - vrl3 = R*(relrot1*ny - relrot2*nx); //- 0.5*((radj-radi)/radsum)*vtr3; - vrlmag = sqrt(vrl1*vrl1+vrl2*vrl2+vrl3*vrl3); - if (vrlmag != 0.0) vrlmaginv = 1.0/vrlmag; - else vrlmaginv = 0.0; - - // Rolling displacement - rollmag = sqrt(shear[3]*shear[3] + shear[4]*shear[4] + shear[5]*shear[5]); - rolldotn = shear[3]*nx + shear[4]*ny + shear[5]*nz; - - if (shearupdate) { - if (fabs(rolldotn) < EPSILON) rolldotn = 0; - if (rolldotn > 0){ //Rotate into tangential plane - scalefac = rollmag/(rollmag - rolldotn); - shear[3] -= rolldotn*nx; - shear[4] -= rolldotn*ny; - shear[5] -= rolldotn*nz; - //Also rescale to preserve magnitude - shear[3] *= scalefac; - shear[4] *= scalefac; - shear[5] *= scalefac; - } - shear[3] += vrl1*dt; - shear[4] += vrl2*dt; - shear[5] += vrl3*dt; - } - - k_R = kR[itype][jtype]*4.0*F_C*pow(aovera0,1.5); - if (rollingdamp == INDEP) eta_R = etaR[itype][jtype]; - else if (rollingdamp == BRILLROLL) eta_R = muR[itype][jtype]*fabs(Fne); - fr1 = -k_R*shear[3] - eta_R*vrl1; - fr2 = -k_R*shear[4] - eta_R*vrl2; - fr3 = -k_R*shear[5] - eta_R*vrl3; - - // rescale frictional displacements and forces if needed - Frcrit = muR[itype][jtype] * fabs(Fne + 2*F_C); - - fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3); - if (fr > Frcrit) { - if (rollmag != 0.0) { - shear[3] = -1.0/k_R*(Frcrit*fr1/fr + eta_R*vrl1); - shear[4] = -1.0/k_R*(Frcrit*fr2/fr + eta_R*vrl2); - shear[5] = -1.0/k_R*(Frcrit*fr3/fr + eta_R*vrl3); - fr1 *= Frcrit/fr; - fr2 *= Frcrit/fr; - fr3 *= Frcrit/fr; - } else fr1 = fr2 = fr3 = 0.0; - } - - - //**************************************** - // Twisting torque, including shear history effects - //**************************************** - magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of Marshall) - shear[6] += magtwist*dt; - k_Q = 0.5*kt*a*a;; //eq 32 - eta_Q = 0.5*eta_T*a*a; - magtortwist = -k_Q*shear[6] - eta_Q*magtwist;//M_t torque (eq 30) - - signtwist = (magtwist > 0) - (magtwist < 0); - Mtcrit=TWOTHIRDS*a*Fscrit;//critical torque (eq 44) - if (fabs(magtortwist) > Mtcrit) { - //shear[6] = Mtcrit/k_Q*magtwist/fabs(magtwist); - shear[6] = 1.0/k_Q*(Mtcrit*signtwist - eta_Q*magtwist); - magtortwist = -Mtcrit * signtwist; //eq 34 - } - - // Apply forces & torques - - fx = nx*Fntot + fs1; - fy = ny*Fntot + fs2; - fz = nz*Fntot + fs3; - - f[i][0] += fx; - f[i][1] += fy; - f[i][2] += fz; - - tor1 = ny*fs3 - nz*fs2; - tor2 = nz*fs1 - nx*fs3; - tor3 = nx*fs2 - ny*fs1; - - torque[i][0] -= radi*tor1; - torque[i][1] -= radi*tor2; - torque[i][2] -= radi*tor3; - - tortwist1 = magtortwist * nx; - tortwist2 = magtortwist * ny; - tortwist3 = magtortwist * nz; - - torque[i][0] += tortwist1; - torque[i][1] += tortwist2; - torque[i][2] += tortwist3; - - torroll1 = R*(ny*fr3 - nz*fr2); //n cross fr - torroll2 = R*(nz*fr1 - nx*fr3); - torroll3 = R*(nx*fr2 - ny*fr1); - - torque[i][0] += torroll1; - torque[i][1] += torroll2; - torque[i][2] += torroll3; - - if (force->newton_pair || j < nlocal) { - f[j][0] -= fx; - f[j][1] -= fy; - f[j][2] -= fz; - - torque[j][0] -= radj*tor1; - torque[j][1] -= radj*tor2; - torque[j][2] -= radj*tor3; - - torque[j][0] -= tortwist1; - torque[j][1] -= tortwist2; - torque[j][2] -= tortwist3; - - torque[j][0] -= torroll1; - torque[j][1] -= torroll2; - torque[j][2] -= torroll3; - } - if (evflag) ev_tally_xyz(i,j,nlocal,0, - 0.0,0.0,fx,fy,fz,delx,dely,delz); - } - } - } -} - -/* ---------------------------------------------------------------------- - global settings - ------------------------------------------------------------------------- */ - -void PairGranJKRRolling::settings(int narg, char **arg) -{ - if (narg < 6) error->all(FLERR,"Illegal pair_style command"); - - int ntypes = atom->ntypes; - - if (narg < 8*ntypes) error->all(FLERR,"Illegal pair_style command"); - - E_one = new double[ntypes+1]; - G_one = new double[ntypes+1]; - pois = new double[ntypes+1]; - muS_one = new double[ntypes+1]; - cor = new double[ntypes+1]; - alpha_one = new double[ntypes+1]; - Ecoh_one = new double[ntypes+1]; - kR_one = new double[ntypes+1]; - muR_one = new double[ntypes+1]; - etaR_one = new double[ntypes+1]; - - //Defaults - normaldamp = TSUJI; - rollingdamp = INDEP; - - int iarg = 8*ntypes; - while (iarg < narg){ - if (strcmp(arg[iarg],"normaldamp") == 0){ - if (iarg+2 > narg) error->all(FLERR, "Invalid pair/gran/dmt/rolling entry"); - if (strcmp(arg[iarg+1],"tsuji") == 0) normaldamp = TSUJI; - else if (strcmp(arg[iarg+1],"brilliantov") == 0) normaldamp = BRILLIANTOV; - else error->all(FLERR, "Invalid normal damping model for pair/gran/dmt/rolling"); - iarg += 2; - } - else if (strcmp(arg[iarg],"rollingdamp") == 0){ - if (iarg+2 > narg) error->all(FLERR, "Invalid pair/gran/dmt/rolling entry"); - if (strcmp(arg[iarg+1],"independent") == 0) rollingdamp = INDEP; - else if (strcmp(arg[iarg+1],"brilliantov") == 0) rollingdamp = BRILLROLL; - else error->all(FLERR, "Invalid rolling damping model for pair/gran/dmt/rolling"); - iarg +=2; - } - else iarg += 1; - } - - for (int i=0; i < ntypes;i++){ - - E_one[i+1] = force->numeric(FLERR, arg[i]); - G_one[i+1] = force->numeric(FLERR, arg[ntypes+i]); - muS_one[i+1] = force->numeric(FLERR, arg[2*ntypes+i]); - cor[i+1] = force->numeric(FLERR, arg[3*ntypes+i]); - Ecoh_one[i+1] = force->numeric(FLERR, arg[4*ntypes+i]); - kR_one[i+1] = force->numeric(FLERR, arg[5*ntypes+i]); - muR_one[i+1] = force->numeric(FLERR, arg[6*ntypes+i]); - etaR_one[i+1] = force->numeric(FLERR, arg[7*ntypes+i]); - } - - //Optional keywords: - // normaldamp tsuji, or normaldamp brilliantov - // rollingdamp brilliantov - - //Derived from inputs - for (int i=1; i <= ntypes; i++){ - pois[i] = E_one[i]/(2.0*G_one[i]) - 1.0; - alpha_one[i] = 1.2728-4.2783*cor[i]+11.087*cor[i]*cor[i]-22.348*cor[i]*cor[i]*cor[i]+27.467*cor[i]*cor[i]*cor[i]*cor[i]-18.022*cor[i]*cor[i]*cor[i]*cor[i]*cor[i]+4.8218*cor[i]*cor[i]*cor[i]*cor[i]*cor[i]*cor[i]; - for (int j=i; j <= ntypes; j++){ - E[i][j] = E[j][i] = 1/((1-pois[i]*pois[i])/E_one[i]+(1-pois[j]*pois[j])/E_one[j]); - G[i][j] = G[j][i] = 1/((2-pois[i])/G_one[i]+(2-pois[j])/G_one[j]); - if (normaldamp == TSUJI){ - alpha[i][j] = alpha[j][i] = sqrt(alpha_one[i]*alpha_one[j]); - } - else if (normaldamp == BRILLIANTOV){ - gamman[i][j] = gamman[j][i] = sqrt(cor[i]*cor[j]); - } - muS[i][j] = muS[j][i] = sqrt(muS_one[i]*muS_one[j]); - Ecoh[i][j] = Ecoh[j][i] = sqrt(Ecoh_one[i]*Ecoh_one[j]); - kR[i][j] = kR[j][i] = sqrt(kR_one[i]*kR_one[j]); - etaR[i][j] = etaR[j][i] = sqrt(etaR_one[i]*etaR_one[j]); - muR[i][j] = muR[j][i] = sqrt(muR_one[i]*muR_one[j]); - } - } -} - -/* ---------------------------------------------------------------------- */ - -double PairGranJKRRolling::single(int i, int j, int itype, int jtype, - double rsq, - double factor_coul, double factor_lj, - double &fforce) -{//TODO: update PairGranJKRRolling::single for JKR - double radi,radj,radsum; - double r,rinv,rsqinv,delx,dely,delz, nx, ny, nz, R; - double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3,wr1,wr2,wr3; - double overlap, a; - double mi,mj,meff,damp,kn,kt; - double Fdamp,Fne,Fntot,Fscrit; - double eta_N,eta_T; - double vtr1,vtr2,vtr3,vrel; - double fs1,fs2,fs3,fs; - double shrmag; - double F_C, delta_C, olapsq, olapcubed, sqrtterm, tmp, a0; - double keyterm, keyterm2, keyterm3, aovera0, foverFc; - - double *radius = atom->radius; - radi = radius[i]; - radj = radius[j]; - radsum = radi + radj; - - r = sqrt(rsq); - rinv = 1.0/r; - rsqinv = 1.0/rsq; - R = radi*radj/(radi+radj); - a0 = pow(9.0*M_PI*Ecoh[itype][jtype]*R*R/E[itype][jtype],ONETHIRD); - delta_C = 0.5*a0*a0*POW6ONE/R; - - int *touch = fix_history->firstflag[i]; - if ((rsq >= (radsum+delta_C)*(radsum+delta_C) )|| - (rsq >= radsum*radsum && touch[j])){ - fforce = 0.0; - svector[0] = svector[1] = svector[2] = svector[3] = 0.0; - return 0.0; - } - - // relative translational velocity - - double **v = atom->v; - vr1 = v[i][0] - v[j][0]; - vr2 = v[i][1] - v[j][1]; - vr3 = v[i][2] - v[j][2]; - - // normal component - - double **x = atom->x; - delx = x[i][0] - x[j][0]; - dely = x[i][1] - x[j][1]; - delz = x[i][2] - x[j][2]; - - nx = delx*rinv; - ny = dely*rinv; - nz = delz*rinv; - - - vnnr = vr1*nx + vr2*ny + vr3*nz; - vn1 = nx*vnnr; - vn2 = ny*vnnr; - vn3 = nz*vnnr; - - // tangential component - - vt1 = vr1 - vn1; - vt2 = vr2 - vn2; - vt3 = vr3 - vn3; - - // relative rotational velocity - - double **omega = atom->omega; - wr1 = (radi*omega[i][0] + radj*omega[j][0]); - wr2 = (radi*omega[i][1] + radj*omega[j][1]); - wr3 = (radi*omega[i][2] + radj*omega[j][2]); - - // meff = effective mass of pair of particles - // if I or J part of rigid body, use body mass - // if I or J is frozen, meff is other particle - - double *rmass = atom->rmass; - int *type = atom->type; - int *mask = atom->mask; - - mi = rmass[i]; - mj = rmass[j]; - if (fix_rigid) { - // NOTE: ensure mass_rigid is current for owned+ghost atoms? - if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; - if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; - } - - meff = mi*mj / (mi+mj); - if (mask[i] & freeze_group_bit) meff = mj; - if (mask[j] & freeze_group_bit) meff = mi; - - - // normal force = JKR - F_C = 3.0*R*M_PI*Ecoh[itype][jtype]; - overlap = radsum - r; - olapsq = overlap*overlap; - olapcubed = olapsq*olapsq; - sqrtterm = sqrt(1.0 + olapcubed); - tmp = 2.0 + olapcubed + 2.0*sqrtterm; - keyterm = pow(tmp,ONETHIRD); - keyterm2 = olapsq/keyterm; - keyterm3 = sqrt(overlap + keyterm2 + keyterm); - aovera0 = POW6TWO * (keyterm3 + - sqrt(2.0*overlap - keyterm2 - keyterm + 4.0/keyterm3));// eq 41 - a = aovera0*a0; - foverFc = 4.0*((aovera0*aovera0*aovera0) - pow(aovera0,1.5));//F_ne/F_C (eq 40) - - Fne = F_C*foverFc; - - //Damping - kn = 4.0/3.0*E[itype][jtype]*a; - if (normaldamp == BRILLIANTOV) eta_N = a*meff*gamman[itype][jtype]; - else if (normaldamp == TSUJI) eta_N=alpha[itype][jtype]*sqrt(meff*kn); - - Fdamp = -eta_N*vnnr; //F_nd eq 23 and Zhao eq 19 - - Fntot = Fne + Fdamp; - - // relative velocities - - vtr1 = vt1 - (nz*wr2-ny*wr3); - vtr2 = vt2 - (nx*wr3-nz*wr1); - vtr3 = vt3 - (ny*wr1-nx*wr2); - vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; - vrel = sqrt(vrel); - - // shear history effects - // neighprev = index of found neigh on previous call - // search entire jnum list of neighbors of I for neighbor J - // start from neighprev, since will typically be next neighbor - // reset neighprev to 0 as necessary - - int jnum = list->numneigh[i]; - int *jlist = list->firstneigh[i]; - double *allshear = fix_history->firstvalue[i]; - - for (int jj = 0; jj < jnum; jj++) { - neighprev++; - if (neighprev >= jnum) neighprev = 0; - if (jlist[neighprev] == j) break; - } - - double *shear = &allshear[size_history*neighprev]; - shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + - shear[2]*shear[2]); - - // tangential forces = shear + tangential velocity damping - kt=8.0*G[itype][jtype]*a; - - eta_T = eta_N; - fs1 = -kt*shear[0] - eta_T*vtr1; - fs2 = -kt*shear[1] - eta_T*vtr2; - fs3 = -kt*shear[2] - eta_T*vtr3; - - // rescale frictional displacements and forces if needed - - fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); - Fscrit= muS[itype][jtype] * fabs(Fne + 2*F_C); - - if (fs > Fscrit) { - if (shrmag != 0.0) { - fs1 *= Fscrit/fs; - fs2 *= Fscrit/fs; - fs3 *= Fscrit/fs; - fs *= Fscrit/fs; - } else fs1 = fs2 = fs3 = fs = 0.0; - } - - // set all forces and return no energy - - fforce = Fntot; - - // set single_extra quantities - - svector[0] = fs1; - svector[1] = fs2; - svector[2] = fs3; - svector[3] = fs; - svector[4] = vn1; - svector[5] = vn2; - svector[6] = vn3; - svector[7] = vt1; - svector[8] = vt2; - svector[9] = vt3; - return 0.0; -} diff --git a/src/GRANULAR/pair_gran_jkr_rolling.h b/src/GRANULAR/pair_gran_jkr_rolling.h deleted file mode 100644 index 8c4b339eb3..0000000000 --- a/src/GRANULAR/pair_gran_jkr_rolling.h +++ /dev/null @@ -1,56 +0,0 @@ -/* -*- c++ -*- ---------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifdef PAIR_CLASS - -PairStyle(gran/jkr/rolling,PairGranJKRRolling) - -#else - -#ifndef LMP_PAIR_GRAN_JKR_ROLLING_H -#define LMP_PAIR_GRAN_JKR_ROLLING_H - -#include "pair_gran_hooke_history.h" - -namespace LAMMPS_NS { - -class PairGranJKRRolling : public PairGranHookeHistory { -public: - PairGranJKRRolling(class LAMMPS *); - virtual ~PairGranJKRRolling(); - virtual void compute(int, int); - void settings(int, char **); //Eventually set this through coeff method so that user can specify a particular i-j set of coefficients - double single(int, int, int, int, double, double, double, double &); - double *E_one, *G_one, *pois, *muS_one, *cor, *alpha_one, *Ecoh_one, *kR_one, *muR_one, *etaR_one; //Public so as to be accessible to fix/wall/gran -private: - double **E, **G, **alpha, **muS, **Ecoh, **kR, **muR, **etaR, **gamman; - int normaldamp, rollingdamp; - - - -}; - -} - -#endif -#endif - -/* ERROR/WARNING messages: - -E: Illegal ... command - -Self-explanatory. Check the input script syntax and compare to the -documentation for the command. You can use -echo screen as a -command-line option when running LAMMPS to see the offending line. - - */ diff --git a/src/GRANULAR/pair_gran_jkr_rolling_multi.cpp b/src/GRANULAR/pair_gran_jkr_rolling_multi.cpp deleted file mode 100644 index a9156390e5..0000000000 --- a/src/GRANULAR/pair_gran_jkr_rolling_multi.cpp +++ /dev/null @@ -1,1181 +0,0 @@ -/* ---------------------------------------------------------------------- -http://lammps.sandia.gov, Sandia National Laboratories -Steve Plimpton, sjplimp@sandia.gov - -Copyright (2003) Sandia Corporation. Under the terms of Contract -DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains -certain rights in this software. This software is distributed under -the GNU General Public License. - -See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - Contributing authors: Leo Silbert (SNL), Gary Grest (SNL) - ------------------------------------------------------------------------- */ - -#include -#include -#include -#include -#include "pair_gran_jkr_rolling_multi.h" -#include "atom.h" -#include "atom_vec.h" -#include "domain.h" -#include "force.h" -#include "update.h" -#include "modify.h" -#include "fix.h" -#include "fix_neigh_history.h" -#include "comm.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "neigh_request.h" -#include "memory.h" -#include "error.h" -#include "math_const.h" -//#include - -using namespace LAMMPS_NS; -using namespace MathConst; - -#define ONETHIRD 0.33333333333333333 -#define TWOTHIRDS 0.66666666666666666 -#define POW6ONE 0.550321208149104 //6^(-1/3) -#define POW6TWO 0.30285343213869 //6^(-2/3) - -#define EPSILON 1e-10 - -enum {TSUJI, BRILLIANTOV}; -enum {INDEP, BRILLROLL}; - -/* ---------------------------------------------------------------------- */ - -PairGranJKRRollingMulti::PairGranJKRRollingMulti(LAMMPS *lmp) : Pair(lmp) -{ - single_enable = 1; - no_virial_fdotr_compute = 1; - history = 1; - fix_history = NULL; - - single_extra = 10; - svector = new double[10]; - - neighprev = 0; - - nmax = 0; - mass_rigid = NULL; - - // set comm size needed by this Pair if used with fix rigid - - comm_forward = 1; -} - -/* ---------------------------------------------------------------------- */ -PairGranJKRRollingMulti::~PairGranJKRRollingMulti() -{ - delete [] svector; - if (fix_history) modify->delete_fix("NEIGH_HISTORY"); - - if (allocated) { - memory->destroy(setflag); - memory->destroy(cutsq); - - memory->destroy(cut); - memory->destroy(E); - memory->destroy(G); - memory->destroy(normaldamp); - memory->destroy(rollingdamp); - memory->destroy(alpha); - memory->destroy(gamman); - memory->destroy(muS); - memory->destroy(Ecoh); - memory->destroy(kR); - memory->destroy(muR); - memory->destroy(etaR); - - delete [] onerad_dynamic; - delete [] onerad_frozen; - delete [] maxrad_dynamic; - delete [] maxrad_frozen; - } - memory->destroy(mass_rigid); -} -/* ---------------------------------------------------------------------- */ - -void PairGranJKRRollingMulti::compute(int eflag, int vflag) -{ -// feenableexcept(FE_INVALID | FE_OVERFLOW); - int i,j,ii,jj,inum,jnum,itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,nx,ny,nz; - double radi,radj,radsum,rsq,r,rinv,rsqinv,R,a; - double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; - double wr1,wr2,wr3; - double vtr1,vtr2,vtr3,vrel; - double kn, kt, k_Q, k_R, eta_N, eta_T, eta_Q, eta_R; - double Fne, Fdamp, Fntot, Fscrit, Frcrit, F_C, delta_C, delta_Cinv; - double overlap, olapsq, olapcubed, sqrtterm, tmp, a0; - double keyterm, keyterm2, keyterm3, aovera0, foverFc; - double mi,mj,meff,damp,ccel,tor1,tor2,tor3; - double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3,vrlmag,vrlmaginv; - double rollmag, rolldotn, scalefac; - double fr, fr1, fr2, fr3; - double signtwist, magtwist, magtortwist, Mtcrit; - double fs,fs1,fs2,fs3,roll1,roll2,roll3,torroll1,torroll2,torroll3; - double tortwist1, tortwist2, tortwist3; - double shrmag,rsht; - int *ilist,*jlist,*numneigh,**firstneigh; - int *touch,**firsttouch; - double *shear,*allshear,**firstshear; - - if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = vflag_fdotr = 0; - - int shearupdate = 1; - if (update->setupflag) shearupdate = 0; - - // update rigid body info for owned & ghost atoms if using FixRigid masses - // body[i] = which body atom I is in, -1 if none - // mass_body = mass of each rigid body - - if (fix_rigid && neighbor->ago == 0){ - int tmp; - int *body = (int *) fix_rigid->extract("body",tmp); - double *mass_body = (double *) fix_rigid->extract("masstotal",tmp); - if (atom->nmax > nmax) { - memory->destroy(mass_rigid); - nmax = atom->nmax; - memory->create(mass_rigid,nmax,"pair:mass_rigid"); - } - int nlocal = atom->nlocal; - for (i = 0; i < nlocal; i++) - if (body[i] >= 0) mass_rigid[i] = mass_body[body[i]]; - else mass_rigid[i] = 0.0; - comm->forward_comm_pair(this); - } - - double **x = atom->x; - double **v = atom->v; - double **f = atom->f; - int *type = atom->type; - double **omega = atom->omega; - double **torque = atom->torque; - double *radius = atom->radius; - double *rmass = atom->rmass; - int *mask = atom->mask; - int nlocal = atom->nlocal; - int newton_pair = force->newton_pair; - - inum = list->inum; - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->firstneigh; - firsttouch = fix_history->firstflag; - firstshear = fix_history->firstvalue; - - // loop over neighbors of my atoms - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - itype = type[i]; - radi = radius[i]; - touch = firsttouch[i]; - allshear = firstshear[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - j &= NEIGHMASK; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - jtype = type[j]; - rsq = delx*delx + dely*dely + delz*delz; - radj = radius[j]; - radsum = radi + radj; - R = radi*radj/(radi+radj); - a0 = pow(9.0*M_PI*Ecoh[itype][jtype]*R*R/E[itype][jtype],ONETHIRD); - delta_C = 0.5*a0*a0*POW6ONE/R; - - if ((rsq >= radsum*radsum && touch[jj] == 0) || - (rsq >= (radsum+delta_C)*(radsum+delta_C))) { - - // unset non-touching neighbors - - touch[jj] = 0; - shear = &allshear[3*jj]; - shear[0] = 0.0; - shear[1] = 0.0; - shear[2] = 0.0; - - } else { - F_C = 3.0*R*M_PI*Ecoh[itype][jtype]; - r = sqrt(rsq); - rinv = 1.0/r; - rsqinv = 1.0/rsq; - - nx = delx*rinv; - ny = dely*rinv; - nz = delz*rinv; - - // relative translational velocity - - vr1 = v[i][0] - v[j][0]; - vr2 = v[i][1] - v[j][1]; - vr3 = v[i][2] - v[j][2]; - - // normal component - - vnnr = vr1*nx + vr2*ny + vr3*nz; //v_R . n - vn1 = nx*vnnr; - vn2 = ny*vnnr; - vn3 = nz*vnnr; - - // meff = effective mass of pair of particles - // if I or J part of rigid body, use body mass - // if I or J is frozen, meff is other particle - - mi = rmass[i]; - mj = rmass[j]; - if (fix_rigid) { - if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; - if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; - } - - meff = mi*mj / (mi+mj); - if (mask[i] & freeze_group_bit) meff = mj; - if (mask[j] & freeze_group_bit) meff = mi; - - //**************************************** - //Normal force = JKR-adjusted Hertzian contact + damping - //**************************************** - if (Ecoh[itype][jtype] != 0.0) delta_Cinv = 1.0/delta_C; - else delta_Cinv = 1.0; - overlap = (radsum - r)*delta_Cinv; - olapsq = overlap*overlap; - olapcubed = olapsq*overlap; - sqrtterm = sqrt(1.0 + olapcubed); - tmp = 2.0 + olapcubed + 2.0*sqrtterm; - keyterm = pow(tmp,ONETHIRD); - keyterm2 = olapsq/keyterm; - keyterm3 = sqrt(overlap + keyterm2 + keyterm); - aovera0 = POW6TWO * (keyterm3 + - sqrt(2.0*overlap - keyterm2 - keyterm + 4.0/keyterm3));// eq 41 - a = aovera0*a0; - foverFc = 4.0*((aovera0*aovera0*aovera0) - pow(aovera0,1.5));//F_ne/F_C (eq 40) - - Fne = F_C*foverFc; - - //Damping - kn = 4.0/3.0*E[itype][jtype]*a; - if (normaldamp[itype][jtype] == BRILLIANTOV) eta_N = a*meff*gamman[itype][jtype]; - else if (normaldamp[itype][jtype] == TSUJI) eta_N=alpha[itype][jtype]*sqrt(meff*kn); - - Fdamp = -eta_N*vnnr; //F_nd eq 23 and Zhao eq 19 - - Fntot = Fne + Fdamp; - //if (screen) fprintf(screen,"%d %d %16.16g %16.16g \n",itype,jtype,Ecoh[itype][jtype],E[itype][jtype]); - //if (logfile) fprintf(logfile,"%d %d %16.16g %16.16g \n",itype,jtype,Ecoh[itype][jtype],E[itype][jtype]); - - //**************************************** - //Tangential force, including shear history effects - //**************************************** - - // tangential component - vt1 = vr1 - vn1; - vt2 = vr2 - vn2; - vt3 = vr3 - vn3; - - // relative rotational velocity - // Luding Gran Matt 2008, v10,p235 suggests correcting radi and radj by subtracting - // delta/2, i.e. instead of radi, use distance to center of contact point? - wr1 = (radi*omega[i][0] + radj*omega[j][0]); - wr2 = (radi*omega[i][1] + radj*omega[j][1]); - wr3 = (radi*omega[i][2] + radj*omega[j][2]); - - // relative tangential velocities - vtr1 = vt1 - (nz*wr2-ny*wr3); - vtr2 = vt2 - (nx*wr3-nz*wr1); - vtr3 = vt3 - (ny*wr1-nx*wr2); - vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; - vrel = sqrt(vrel); - - // shear history effects - touch[jj] = 1; - shear = &allshear[3*jj]; - shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + - shear[2]*shear[2]); - - // Rotate and update shear displacements. - // See e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 - if (shearupdate) { - rsht = shear[0]*nx + shear[1]*ny + shear[2]*nz; - if (fabs(rsht) < EPSILON) rsht = 0; - if (rsht > 0){ - scalefac = shrmag/(shrmag - rsht); //if rhst == shrmag, contacting pair has rotated 90 deg. in one step, in which case you deserve a crash! - shear[0] -= rsht*nx; - shear[1] -= rsht*ny; - shear[2] -= rsht*nz; - //Also rescale to preserve magnitude - shear[0] *= scalefac; - shear[1] *= scalefac; - shear[2] *= scalefac; - } - //Update shear history - shear[0] += vtr1*dt; - shear[1] += vtr2*dt; - shear[2] += vtr3*dt; - } - - // tangential forces = shear + tangential velocity damping - // following Zhao and Marshall Phys Fluids v20, p043302 (2008) - kt=8.0*G[itype][jtype]*a; - - eta_T = eta_N; //Based on discussion in Marshall; eta_T can also be an independent parameter - fs1 = -kt*shear[0] - eta_T*vtr1; //eq 26 - fs2 = -kt*shear[1] - eta_T*vtr2; - fs3 = -kt*shear[2] - eta_T*vtr3; - - // rescale frictional displacements and forces if needed - Fscrit = muS[itype][jtype] * fabs(Fne + 2*F_C); - // For JKR, use eq 43 of Marshall. For DMT, use Fne instead - - fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); - if (fs > Fscrit) { - if (shrmag != 0.0) { - //shear[0] = (Fcrit/fs) * (shear[0] + eta_T*vtr1/kt) - eta_T*vtr1/kt; - //shear[1] = (Fcrit/fs) * (shear[1] + eta_T*vtr1/kt) - eta_T*vtr1/kt; - //shear[2] = (Fcrit/fs) * (shear[2] + eta_T*vtr1/kt) - eta_T*vtr1/kt; - shear[0] = -1.0/kt*(Fscrit*fs1/fs + eta_T*vtr1); //Same as above, but simpler (check!) - shear[1] = -1.0/kt*(Fscrit*fs2/fs + eta_T*vtr2); - shear[2] = -1.0/kt*(Fscrit*fs3/fs + eta_T*vtr3); - fs1 *= Fscrit/fs; - fs2 *= Fscrit/fs; - fs3 *= Fscrit/fs; - } else fs1 = fs2 = fs3 = 0.0; - } - - //**************************************** - // Rolling force, including shear history effects - //**************************************** - - relrot1 = omega[i][0] - omega[j][0]; - relrot2 = omega[i][1] - omega[j][1]; - relrot3 = omega[i][2] - omega[j][2]; - - // rolling velocity, see eq. 31 of Wang et al, Particuology v 23, p 49 (2015) - // This is different from the Marshall papers, which use the Bagi/Kuhn formulation - // for rolling velocity (see Wang et al for why the latter is wrong) - vrl1 = R*(relrot2*nz - relrot3*ny); //- 0.5*((radj-radi)/radsum)*vtr1; - vrl2 = R*(relrot3*nx - relrot1*nz); //- 0.5*((radj-radi)/radsum)*vtr2; - vrl3 = R*(relrot1*ny - relrot2*nx); //- 0.5*((radj-radi)/radsum)*vtr3; - vrlmag = sqrt(vrl1*vrl1+vrl2*vrl2+vrl3*vrl3); - if (vrlmag != 0.0) vrlmaginv = 1.0/vrlmag; - else vrlmaginv = 0.0; - - // Rolling displacement - rollmag = sqrt(shear[3]*shear[3] + shear[4]*shear[4] + shear[5]*shear[5]); - rolldotn = shear[3]*nx + shear[4]*ny + shear[5]*nz; - - if (shearupdate) { - if (fabs(rolldotn) < EPSILON) rolldotn = 0; - if (rolldotn > 0){ //Rotate into tangential plane - scalefac = rollmag/(rollmag - rolldotn); - shear[3] -= rolldotn*nx; - shear[4] -= rolldotn*ny; - shear[5] -= rolldotn*nz; - //Also rescale to preserve magnitude - shear[3] *= scalefac; - shear[4] *= scalefac; - shear[5] *= scalefac; - } - shear[3] += vrl1*dt; - shear[4] += vrl2*dt; - shear[5] += vrl3*dt; - } - - k_R = kR[itype][jtype]*4.0*F_C*pow(aovera0,1.5); - if (rollingdamp[itype][jtype] == INDEP) eta_R = etaR[itype][jtype]; - else if (rollingdamp[itype][jtype] == BRILLROLL) eta_R = muR[itype][jtype]*fabs(Fne); - fr1 = -k_R*shear[3] - eta_R*vrl1; - fr2 = -k_R*shear[4] - eta_R*vrl2; - fr3 = -k_R*shear[5] - eta_R*vrl3; - - // rescale frictional displacements and forces if needed - Frcrit = muR[itype][jtype] * fabs(Fne + 2*F_C); - - fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3); - if (fr > Frcrit) { - if (rollmag != 0.0) { - shear[3] = -1.0/k_R*(Frcrit*fr1/fr + eta_R*vrl1); - shear[4] = -1.0/k_R*(Frcrit*fr2/fr + eta_R*vrl2); - shear[5] = -1.0/k_R*(Frcrit*fr3/fr + eta_R*vrl3); - fr1 *= Frcrit/fr; - fr2 *= Frcrit/fr; - fr3 *= Frcrit/fr; - } else fr1 = fr2 = fr3 = 0.0; - } - - - //**************************************** - // Twisting torque, including shear history effects - //**************************************** - magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of Marshall) - shear[6] += magtwist*dt; - k_Q = 0.5*kt*a*a;; //eq 32 - eta_Q = 0.5*eta_T*a*a; - magtortwist = -k_Q*shear[6] - eta_Q*magtwist;//M_t torque (eq 30) - - signtwist = (magtwist > 0) - (magtwist < 0); - Mtcrit=TWOTHIRDS*a*Fscrit;//critical torque (eq 44) - if (fabs(magtortwist) > Mtcrit) { - //shear[6] = Mtcrit/k_Q*magtwist/fabs(magtwist); - shear[6] = 1.0/k_Q*(Mtcrit*signtwist - eta_Q*magtwist); - magtortwist = -Mtcrit * signtwist; //eq 34 - } - - // Apply forces & torques - - fx = nx*Fntot + fs1; - fy = ny*Fntot + fs2; - fz = nz*Fntot + fs3; - - //if (screen) fprintf(screen,"%16.16g %16.16g %16.16g %16.16g %16.16g %16.16g %16.16g \n",fs1,fs2,fs3,Fntot,nx,ny,nz); - //if (logfile) fprintf(logfile,"%16.16g %16.16g %16.16g %16.16g %16.16g %16.16g %16.16g \n",fs1,fs2,fs3,Fntot,nx,ny,nz); - - f[i][0] += fx; - f[i][1] += fy; - f[i][2] += fz; - - tor1 = ny*fs3 - nz*fs2; - tor2 = nz*fs1 - nx*fs3; - tor3 = nx*fs2 - ny*fs1; - - torque[i][0] -= radi*tor1; - torque[i][1] -= radi*tor2; - torque[i][2] -= radi*tor3; - - tortwist1 = magtortwist * nx; - tortwist2 = magtortwist * ny; - tortwist3 = magtortwist * nz; - - torque[i][0] += tortwist1; - torque[i][1] += tortwist2; - torque[i][2] += tortwist3; - - torroll1 = R*(ny*fr3 - nz*fr2); //n cross fr - torroll2 = R*(nz*fr1 - nx*fr3); - torroll3 = R*(nx*fr2 - ny*fr1); - - torque[i][0] += torroll1; - torque[i][1] += torroll2; - torque[i][2] += torroll3; - - if (force->newton_pair || j < nlocal) { - f[j][0] -= fx; - f[j][1] -= fy; - f[j][2] -= fz; - - torque[j][0] -= radj*tor1; - torque[j][1] -= radj*tor2; - torque[j][2] -= radj*tor3; - - torque[j][0] -= tortwist1; - torque[j][1] -= tortwist2; - torque[j][2] -= tortwist3; - - torque[j][0] -= torroll1; - torque[j][1] -= torroll2; - torque[j][2] -= torroll3; - } - if (evflag) ev_tally_xyz(i,j,nlocal,0, - 0.0,0.0,fx,fy,fz,delx,dely,delz); - } - } - } -} - - -/* ---------------------------------------------------------------------- - allocate all arrays - ------------------------------------------------------------------------- */ - -void PairGranJKRRollingMulti::allocate() -{ - allocated = 1; - int n = atom->ntypes; - - memory->create(setflag,n+1,n+1,"pair:setflag"); - for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) - setflag[i][j] = 0; - - memory->create(cutsq,n+1,n+1,"pair:cutsq"); - memory->create(cut,n+1,n+1,"pair:cut"); - memory->create(E,n+1,n+1,"pair:E"); - memory->create(G,n+1,n+1,"pair:G"); - memory->create(normaldamp,n+1,n+1,"pair:normaldamp"); - memory->create(rollingdamp,n+1,n+1,"pair:rollingdamp"); - memory->create(alpha,n+1,n+1,"pair:alpha"); - memory->create(gamman,n+1,n+1,"pair:gamman"); - memory->create(muS,n+1,n+1,"pair:muS"); - memory->create(Ecoh,n+1,n+1,"pair:Ecoh"); - memory->create(kR,n+1,n+1,"pair:kR"); - memory->create(muR,n+1,n+1,"pair:muR"); - memory->create(etaR,n+1,n+1,"pair:etaR"); - - onerad_dynamic = new double[n+1]; - onerad_frozen = new double[n+1]; - maxrad_dynamic = new double[n+1]; - maxrad_frozen = new double[n+1]; -} - -/* ---------------------------------------------------------------------- - global settings - ------------------------------------------------------------------------- */ - -void PairGranJKRRollingMulti::settings(int narg, char **arg) -{ - if (narg != 1) error->all(FLERR,"Illegal pair_style command"); - - if (strcmp(arg[0],"NULL") == 0 ) cut_global = -1.0; - else cut_global = force->numeric(FLERR,arg[0]); - - // reset cutoffs that have been explicitly set - if (allocated) { - int i,j; - for (i = 1; i <= atom->ntypes; i++) - for (j = i; j <= atom->ntypes; j++) - if (setflag[i][j]) cut[i][j] = cut_global; - } -} - -/* ---------------------------------------------------------------------- - set coeffs for one or more type pairs - ------------------------------------------------------------------------- */ - -void PairGranJKRRollingMulti::coeff(int narg, char **arg) -{ - if (narg < 10 || narg > 15) - error->all(FLERR,"Incorrect args for pair coefficients2"); - - if (!allocated) allocate(); - - int ilo,ihi,jlo,jhi; - force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi); - force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi); - - double E_one = force->numeric(FLERR,arg[2]); - double G_one = force->numeric(FLERR,arg[3]); - double muS_one = force->numeric(FLERR,arg[4]); - double cor_one = force->numeric(FLERR,arg[5]); - double Ecoh_one = force->numeric(FLERR,arg[6]); - double kR_one = force->numeric(FLERR,arg[7]); - double muR_one = force->numeric(FLERR,arg[8]); - double etaR_one = force->numeric(FLERR,arg[9]); - - //Defaults - int normaldamp_one = TSUJI; - int rollingdamp_one = INDEP; - double cut_one = cut_global; - - int iarg = 10; - while (iarg < narg) { - if (strcmp(arg[iarg],"normaldamp") == 0){ - if (iarg+2 > narg) error->all(FLERR, "Invalid pair/gran/dmt/rolling entry"); - if (strcmp(arg[iarg+1],"tsuji") == 0) normaldamp_one = TSUJI; - else if (strcmp(arg[iarg+1],"brilliantov") == 0) normaldamp_one = BRILLIANTOV; - else error->all(FLERR, "Invalid normal damping model for pair/gran/dmt/rolling"); - iarg += 2; - } - else if (strcmp(arg[iarg],"rollingdamp") == 0){ - if (iarg+2 > narg) error->all(FLERR, "Invalid pair/gran/dmt/rolling entry"); - if (strcmp(arg[iarg+1],"independent") == 0) rollingdamp_one = INDEP; - else if (strcmp(arg[iarg+1],"brilliantov") == 0) rollingdamp_one = BRILLROLL; - else error->all(FLERR, "Invalid rolling damping model for pair/gran/dmt/rolling"); - iarg +=2; - } - else { - if (strcmp(arg[iarg],"NULL") == 0) cut_one = -1.0; - else cut_one = force->numeric(FLERR,arg[iarg]); - iarg += 1; - } - } - - int count = 0; - for (int i = ilo; i <= ihi; i++) { - double pois = E_one/(2.0*G_one) - 1.0; - double alpha_one = 1.2728-4.2783*cor_one+11.087*cor_one*cor_one-22.348*cor_one*cor_one*cor_one+27.467*cor_one*cor_one*cor_one*cor_one-18.022*cor_one*cor_one*cor_one*cor_one*cor_one+4.8218*cor_one*cor_one*cor_one*cor_one*cor_one*cor_one; - - for (int j = MAX(jlo,i); j <= jhi; j++) { - E[i][j] = E_one; - G[i][j] = G_one; - if (normaldamp_one == TSUJI) { - normaldamp[i][j] = TSUJI; - alpha[i][j] = alpha_one; - } - else if (normaldamp_one == BRILLIANTOV) { - normaldamp[i][j] = BRILLIANTOV; - gamman[i][j] = cor_one; - } - if (rollingdamp_one == INDEP) { - rollingdamp[i][j] = INDEP; - } - else if (rollingdamp_one == BRILLROLL) { - rollingdamp[i][j] = BRILLROLL; - } - muS[i][j] = muS_one; - Ecoh[i][j] = Ecoh_one; - kR[i][j] = kR_one; - etaR[i][j] = etaR_one; - muR[i][j] = muR_one; - cut[i][j] = cut_one; - setflag[i][j] = 1; - count++; - } - } - - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients1"); -} - -/* ---------------------------------------------------------------------- - init specific to this pair style - ------------------------------------------------------------------------- */ - -void PairGranJKRRollingMulti::init_style() -{ - int i; - - // error and warning checks - - if (!atom->radius_flag || !atom->rmass_flag) - error->all(FLERR,"Pair granular requires atom attributes radius, rmass"); - if (comm->ghost_velocity == 0) - error->all(FLERR,"Pair granular requires ghost atoms store velocity"); - - // need a granular neigh list - - int irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->size = 1; - if (history) neighbor->requests[irequest]->history = 1; - - dt = update->dt; - - // if shear history is stored: - // if first init, create Fix needed for storing shear history - - if (history && fix_history == NULL) { - char dnumstr[16]; - sprintf(dnumstr,"%d",3); - char **fixarg = new char*[4]; - fixarg[0] = (char *) "NEIGH_HISTORY"; - fixarg[1] = (char *) "all"; - fixarg[2] = (char *) "NEIGH_HISTORY"; - fixarg[3] = dnumstr; - modify->add_fix(4,fixarg,1); - delete [] fixarg; - fix_history = (FixNeighHistory *) modify->fix[modify->nfix-1]; - fix_history->pair = this; - } - - // check for FixFreeze and set freeze_group_bit - - for (i = 0; i < modify->nfix; i++) - if (strcmp(modify->fix[i]->style,"freeze") == 0) break; - if (i < modify->nfix) freeze_group_bit = modify->fix[i]->groupbit; - else freeze_group_bit = 0; - - // check for FixRigid so can extract rigid body masses - - fix_rigid = NULL; - for (i = 0; i < modify->nfix; i++) - if (modify->fix[i]->rigid_flag) break; - if (i < modify->nfix) fix_rigid = modify->fix[i]; - - // check for FixPour and FixDeposit so can extract particle radii - - int ipour; - for (ipour = 0; ipour < modify->nfix; ipour++) - if (strcmp(modify->fix[ipour]->style,"pour") == 0) break; - if (ipour == modify->nfix) ipour = -1; - - int idep; - for (idep = 0; idep < modify->nfix; idep++) - if (strcmp(modify->fix[idep]->style,"deposit") == 0) break; - if (idep == modify->nfix) idep = -1; - - // set maxrad_dynamic and maxrad_frozen for each type - // include future FixPour and FixDeposit particles as dynamic - - int itype; - for (i = 1; i <= atom->ntypes; i++) { - onerad_dynamic[i] = onerad_frozen[i] = 0.0; - if (ipour >= 0) { - itype = i; - onerad_dynamic[i] = - *((double *) modify->fix[ipour]->extract("radius",itype)); - } - if (idep >= 0) { - itype = i; - onerad_dynamic[i] = - *((double *) modify->fix[idep]->extract("radius",itype)); - } - } - - double *radius = atom->radius; - int *mask = atom->mask; - int *type = atom->type; - int nlocal = atom->nlocal; - - for (i = 0; i < nlocal; i++) - if (mask[i] & freeze_group_bit) - onerad_frozen[type[i]] = MAX(onerad_frozen[type[i]],radius[i]); - else - onerad_dynamic[type[i]] = MAX(onerad_dynamic[type[i]],radius[i]); - - MPI_Allreduce(&onerad_dynamic[1],&maxrad_dynamic[1],atom->ntypes, - MPI_DOUBLE,MPI_MAX,world); - MPI_Allreduce(&onerad_frozen[1],&maxrad_frozen[1],atom->ntypes, - MPI_DOUBLE,MPI_MAX,world); - - // set fix which stores history info - - if (history) { - int ifix = modify->find_fix("NEIGH_HISTORY"); - if (ifix < 0) error->all(FLERR,"Could not find pair fix neigh history ID"); - fix_history = (FixNeighHistory *) modify->fix[ifix]; - } -} - -/* ---------------------------------------------------------------------- - init for one type pair i,j and corresponding j,i - ------------------------------------------------------------------------- */ - -double PairGranJKRRollingMulti::init_one(int i, int j) -{ - if (setflag[i][j] == 0) { - E[i][j] = mix_stiffnessE(E[i][i],E[j][j],G[i][i],G[j][j]); - G[i][j] = mix_stiffnessG(G[i][i],E[j][j],G[i][i],G[j][j]); - if (normaldamp[i][j] == TSUJI) { - alpha[i][j] = mix_geom(alpha[i][i],alpha[j][j]); - } - else if (normaldamp[i][j] == BRILLIANTOV) { - gamman[i][j] = mix_geom(gamman[i][i],gamman[j][j]); - } - muS[i][j] = mix_geom(muS[i][i],muS[j][j]); - Ecoh[i][j] = mix_geom(Ecoh[i][i],Ecoh[j][j]); - kR[i][j] = mix_geom(kR[i][i],kR[j][j]); - etaR[i][j] = mix_geom(etaR[i][i],etaR[j][j]); - muR[i][j] = mix_geom(muR[i][i],muR[j][j]); - } - - E[j][i] = E[i][j]; - G[j][i] = G[i][j]; - normaldamp[j][i] = normaldamp[i][j]; - alpha[j][i] = alpha[i][j]; - gamman[j][i] = gamman[i][j]; - rollingdamp[j][i] = rollingdamp[i][j]; - muS[j][i] = muS[i][j]; - Ecoh[j][i] = Ecoh[i][j]; - kR[j][i] = kR[i][j]; - etaR[j][i] = etaR[i][j]; - muR[j][i] = muR[i][j]; - - double cutoff = cut[i][j]; - - // It is likely that cut[i][j] at this point is still 0.0. This can happen when - // there is a future fix_pour after the current run. A cut[i][j] = 0.0 creates - // problems because neighbor.cpp uses min(cut[i][j]) to decide on the bin size - // To avoid this issue,for cases involving cut[i][j] = 0.0 (possible only - // if there is no current information about radius/cutoff of type i and j). - // we assign cutoff = min(cut[i][j]) for i,j such that cut[i][j] > 0.0. - - if (cut[i][j] < 0.0) { - if (((maxrad_dynamic[i] > 0.0) && (maxrad_dynamic[j] > 0.0)) || ((maxrad_dynamic[i] > 0.0) && (maxrad_frozen[j] > 0.0)) || - ((maxrad_frozen[i] > 0.0) && (maxrad_dynamic[j] > 0.0))) { // radius info about both i and j exist - cutoff = maxrad_dynamic[i]+maxrad_dynamic[j]; - cutoff = MAX(cutoff,maxrad_frozen[i]+maxrad_dynamic[j]); - cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]); - } - else { // radius info about either i or j does not exist (i.e. not present and not about to get poured; set to largest value to not interfere with neighbor list) - double cutmax = 0.0; - for (int k = 1; k <= atom->ntypes; k++) { - cutmax = MAX(cutmax,2.0*maxrad_dynamic[k]); - cutmax = MAX(cutmax,2.0*maxrad_frozen[k]); - } - cutoff = cutmax; - } - } - return cutoff; -} - - -/* ---------------------------------------------------------------------- - proc 0 writes to restart file - ------------------------------------------------------------------------- */ - -void PairGranJKRRollingMulti::write_restart(FILE *fp) -{ - write_restart_settings(fp); - - int i,j; - for (i = 1; i <= atom->ntypes; i++) { - for (j = i; j <= atom->ntypes; j++) { - fwrite(&setflag[i][j],sizeof(int),1,fp); - if (setflag[i][j]) { - fwrite(&E[i][j],sizeof(double),1,fp); - fwrite(&G[i][j],sizeof(double),1,fp); - fwrite(&normaldamp[i][j],sizeof(int),1,fp); - fwrite(&rollingdamp[i][j],sizeof(int),1,fp); - fwrite(&alpha[i][j],sizeof(double),1,fp); - fwrite(&gamman[i][j],sizeof(double),1,fp); - fwrite(&muS[i][j],sizeof(double),1,fp); - fwrite(&Ecoh[i][j],sizeof(double),1,fp); - fwrite(&kR[i][j],sizeof(double),1,fp); - fwrite(&muR[i][j],sizeof(double),1,fp); - fwrite(&etaR[i][j],sizeof(double),1,fp); - fwrite(&cut[i][j],sizeof(double),1,fp); - } - } - } -} - -/* ---------------------------------------------------------------------- - proc 0 reads from restart file, bcasts - ------------------------------------------------------------------------- */ - -void PairGranJKRRollingMulti::read_restart(FILE *fp) -{ - read_restart_settings(fp); - allocate(); - - int i,j; - int me = comm->me; - for (i = 1; i <= atom->ntypes; i++) { - for (j = i; j <= atom->ntypes; j++) { - if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); - MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); - if (setflag[i][j]) { - if (me == 0) { - fread(&E[i][j],sizeof(double),1,fp); - fread(&G[i][j],sizeof(double),1,fp); - fread(&normaldamp[i][j],sizeof(int),1,fp); - fread(&rollingdamp[i][j],sizeof(int),1,fp); - fread(&alpha[i][j],sizeof(double),1,fp); - fread(&gamman[i][j],sizeof(double),1,fp); - fread(&muS[i][j],sizeof(double),1,fp); - fread(&Ecoh[i][j],sizeof(double),1,fp); - fread(&kR[i][j],sizeof(double),1,fp); - fread(&muR[i][j],sizeof(double),1,fp); - fread(&etaR[i][j],sizeof(double),1,fp); - fread(&cut[i][j],sizeof(double),1,fp); - } - MPI_Bcast(&E[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&G[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&normaldamp[i][j],1,MPI_INT,0,world); - MPI_Bcast(&rollingdamp[i][j],1,MPI_INT,0,world); - MPI_Bcast(&alpha[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&gamman[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&muS[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&Ecoh[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&kR[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&muR[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&etaR[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); - } - } - } -} - -/* ---------------------------------------------------------------------- - proc 0 writes to restart file - ------------------------------------------------------------------------- */ - -void PairGranJKRRollingMulti::write_restart_settings(FILE *fp) -{ - fwrite(&cut_global,sizeof(double),1,fp); -} - -/* ---------------------------------------------------------------------- - proc 0 reads from restart file, bcasts - ------------------------------------------------------------------------- */ - -void PairGranJKRRollingMulti::read_restart_settings(FILE *fp) -{ - if (comm->me == 0) { - fread(&cut_global,sizeof(double),1,fp); - } - MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); -} - -/* ---------------------------------------------------------------------- */ - -void PairGranJKRRollingMulti::reset_dt() -{ - dt = update->dt; -} - -/* ---------------------------------------------------------------------- */ - -double PairGranJKRRollingMulti::single(int i, int j, int itype, int jtype, - double rsq, double factor_coul, double factor_lj, double &fforce) -{ -// feenableexcept(FE_INVALID | FE_OVERFLOW); - double radi,radj,radsum; - double r,rinv,rsqinv,delx,dely,delz, nx, ny, nz, R; - double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3,wr1,wr2,wr3; - double overlap, a; - double mi,mj,meff,damp,kn,kt; - double Fdamp,Fne,Fntot,Fscrit; - double eta_N,eta_T; - double vtr1,vtr2,vtr3,vrel; - double fs1,fs2,fs3,fs; - double shrmag; - double F_C, delta_C, olapsq, olapcubed, sqrtterm, tmp, a0; - double keyterm, keyterm2, keyterm3, aovera0, foverFc; - - double *radius = atom->radius; - radi = radius[i]; - radj = radius[j]; - radsum = radi + radj; - - r = sqrt(rsq); - rinv = 1.0/r; - rsqinv = 1.0/rsq; - R = radi*radj/(radi+radj); - a0 = pow(9.0*M_PI*Ecoh[itype][jtype]*R*R/E[itype][jtype],ONETHIRD); - delta_C = 0.5*a0*a0*POW6ONE/R; - - int *touch = fix_history->firstflag[i]; - if ((rsq >= (radsum+delta_C)*(radsum+delta_C) )|| - (rsq >= radsum*radsum && touch[j])){ - fforce = 0.0; - svector[0] = svector[1] = svector[2] = svector[3] = 0.0; - return 0.0; - } - - // relative translational velocity - - double **v = atom->v; - vr1 = v[i][0] - v[j][0]; - vr2 = v[i][1] - v[j][1]; - vr3 = v[i][2] - v[j][2]; - - // normal component - - double **x = atom->x; - delx = x[i][0] - x[j][0]; - dely = x[i][1] - x[j][1]; - delz = x[i][2] - x[j][2]; - - nx = delx*rinv; - ny = dely*rinv; - nz = delz*rinv; - - - vnnr = vr1*nx + vr2*ny + vr3*nz; - vn1 = nx*vnnr; - vn2 = ny*vnnr; - vn3 = nz*vnnr; - - // tangential component - - vt1 = vr1 - vn1; - vt2 = vr2 - vn2; - vt3 = vr3 - vn3; - - // relative rotational velocity - - double **omega = atom->omega; - wr1 = (radi*omega[i][0] + radj*omega[j][0]); - wr2 = (radi*omega[i][1] + radj*omega[j][1]); - wr3 = (radi*omega[i][2] + radj*omega[j][2]); - - // meff = effective mass of pair of particles - // if I or J part of rigid body, use body mass - // if I or J is frozen, meff is other particle - - double *rmass = atom->rmass; - int *type = atom->type; - int *mask = atom->mask; - - mi = rmass[i]; - mj = rmass[j]; - if (fix_rigid) { - // NOTE: ensure mass_rigid is current for owned+ghost atoms? - if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; - if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; - } - - meff = mi*mj / (mi+mj); - if (mask[i] & freeze_group_bit) meff = mj; - if (mask[j] & freeze_group_bit) meff = mi; - - - // normal force = JKR - F_C = 3.0*R*M_PI*Ecoh[itype][jtype]; - overlap = radsum - r; - olapsq = overlap*overlap; - olapcubed = olapsq*olapsq; - sqrtterm = sqrt(1.0 + olapcubed); - tmp = 2.0 + olapcubed + 2.0*sqrtterm; - keyterm = pow(tmp,ONETHIRD); - keyterm2 = olapsq/keyterm; - keyterm3 = sqrt(overlap + keyterm2 + keyterm); - aovera0 = POW6TWO * (keyterm3 + - sqrt(2.0*overlap - keyterm2 - keyterm + 4.0/keyterm3));// eq 41 - a = aovera0*a0; - foverFc = 4.0*((aovera0*aovera0*aovera0) - pow(aovera0,1.5));//F_ne/F_C (eq 40) - - Fne = F_C*foverFc; - - //Damping - kn = 4.0/3.0*E[itype][jtype]*a; - if (normaldamp[itype][jtype] == BRILLIANTOV) eta_N = a*meff*gamman[itype][jtype]; - else if (normaldamp[itype][jtype] == TSUJI) eta_N=alpha[itype][jtype]*sqrt(meff*kn); - - Fdamp = -eta_N*vnnr; //F_nd eq 23 and Zhao eq 19 - - Fntot = Fne + Fdamp; - - // relative velocities - - vtr1 = vt1 - (nz*wr2-ny*wr3); - vtr2 = vt2 - (nx*wr3-nz*wr1); - vtr3 = vt3 - (ny*wr1-nx*wr2); - vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; - vrel = sqrt(vrel); - - // shear history effects - // neighprev = index of found neigh on previous call - // search entire jnum list of neighbors of I for neighbor J - // start from neighprev, since will typically be next neighbor - // reset neighprev to 0 as necessary - - int jnum = list->numneigh[i]; - int *jlist = list->firstneigh[i]; - double *allshear = fix_history->firstvalue[i]; - - for (int jj = 0; jj < jnum; jj++) { - neighprev++; - if (neighprev >= jnum) neighprev = 0; - if (jlist[neighprev] == j) break; - } - - double *shear = &allshear[3*neighprev]; - shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + - shear[2]*shear[2]); - - // tangential forces = shear + tangential velocity damping - kt=8.0*G[itype][jtype]*a; - - eta_T = eta_N; - fs1 = -kt*shear[0] - eta_T*vtr1; - fs2 = -kt*shear[1] - eta_T*vtr2; - fs3 = -kt*shear[2] - eta_T*vtr3; - - // rescale frictional displacements and forces if needed - - fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); - Fscrit= muS[itype][jtype] * fabs(Fne + 2*F_C); - - if (fs > Fscrit) { - if (shrmag != 0.0) { - fs1 *= Fscrit/fs; - fs2 *= Fscrit/fs; - fs3 *= Fscrit/fs; - fs *= Fscrit/fs; - } else fs1 = fs2 = fs3 = fs = 0.0; - } - - // set all forces and return no energy - - fforce = Fntot; - - // set single_extra quantities - - svector[0] = fs1; - svector[1] = fs2; - svector[2] = fs3; - svector[3] = fs; - svector[4] = vn1; - svector[5] = vn2; - svector[6] = vn3; - svector[7] = vt1; - svector[8] = vt2; - svector[9] = vt3; - return 0.0; -} - -/* ---------------------------------------------------------------------- */ - -int PairGranJKRRollingMulti::pack_forward_comm(int n, int *list, double *buf, - int pbc_flag, int *pbc) -{ - int i,j,m; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = mass_rigid[j]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -void PairGranJKRRollingMulti::unpack_forward_comm(int n, int first, double *buf) -{ - int i,m,last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) - mass_rigid[i] = buf[m++]; -} - -/* ---------------------------------------------------------------------- - memory usage of local atom-based arrays - ------------------------------------------------------------------------- */ - -double PairGranJKRRollingMulti::memory_usage() -{ - double bytes = nmax * sizeof(double); - return bytes; -} - -/* ---------------------------------------------------------------------- - mixing of stiffness (E) - ------------------------------------------------------------------------- */ - -double PairGranJKRRollingMulti::mix_stiffnessE(double Eii, double Ejj, double Gii, double Gjj) -{ - double poisii = Eii/(2.0*Gii) - 1.0; - double poisjj = Ejj/(2.0*Gjj) - 1.0; - return 1/((1-poisii*poisjj)/Eii+(1-poisjj*poisjj)/Ejj); -} - -/* ---------------------------------------------------------------------- - mixing of stiffness (G) - ------------------------------------------------------------------------- */ - -double PairGranJKRRollingMulti::mix_stiffnessG(double Eii, double Ejj, double Gii, double Gjj) -{ - double poisii = Eii/(2.0*Gii) - 1.0; - double poisjj = Ejj/(2.0*Gjj) - 1.0; - return 1/((2.0 -poisjj)/Gii+(2.0-poisjj)/Gjj); -} - -/* ---------------------------------------------------------------------- - mixing of everything else - ------------------------------------------------------------------------- */ - -double PairGranJKRRollingMulti::mix_geom(double valii, double valjj) -{ - return sqrt(valii*valjj); -} diff --git a/src/GRANULAR/pair_gran_jkr_rolling_multi.h b/src/GRANULAR/pair_gran_jkr_rolling_multi.h deleted file mode 100644 index c9c75de9a6..0000000000 --- a/src/GRANULAR/pair_gran_jkr_rolling_multi.h +++ /dev/null @@ -1,87 +0,0 @@ -/* -*- c++ -*- ---------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifdef PAIR_CLASS - -PairStyle(gran/jkr/rolling/multi,PairGranJKRRollingMulti) - -#else - -#ifndef LMP_PAIR_GRAN_JKR_ROLLING_MULTI_H -#define LMP_PAIR_GRAN_JKR_ROLLING_MULTI_H - -#include "pair.h" - -namespace LAMMPS_NS { - -class PairGranJKRRollingMulti : public Pair { -public: - PairGranJKRRollingMulti(class LAMMPS *); - virtual ~PairGranJKRRollingMulti(); - virtual void compute(int, int); - virtual void settings(int, char **); - virtual void coeff(int, char **); // Made Virtual by IS Oct 7 2017 - void init_style(); - double init_one(int, int); - void write_restart(FILE *); - void read_restart(FILE *); - void write_restart_settings(FILE *); - void read_restart_settings(FILE *); - void reset_dt(); - virtual double single(int, int, int, int, double, double, double, double &); - int pack_forward_comm(int, int *, double *, int, int *); - void unpack_forward_comm(int, int, double *); - double memory_usage(); - - protected: - double cut_global; - double **E,**G,**alpha,**gamman,**muS,**Ecoh,**kR,**muR,**etaR,**cut; - int **normaldamp, **rollingdamp; - double dt; - int freeze_group_bit; - int history; - - int neighprev; - double *onerad_dynamic,*onerad_frozen; - double *maxrad_dynamic,*maxrad_frozen; - - class FixNeighHistory *fix_history; - - // storage of rigid body masses for use in granular interactions - - class Fix *fix_rigid; // ptr to rigid body fix, NULL if none - double *mass_rigid; // rigid mass for owned+ghost atoms - int nmax; // allocated size of mass_rigid - - virtual void allocate(); // Made Virtual by IS Oct 7 2017 - -private: - double mix_stiffnessE(double Eii, double Ejj, double Gii, double Gjj); - double mix_stiffnessG(double Eii, double Ejj, double Gii, double Gjj); - double mix_geom(double valii, double valjj); -}; - -} - -#endif -#endif - -/* ERROR/WARNING messages: - -E: Illegal ... command - -Self-explanatory. Check the input script syntax and compare to the -documentation for the command. You can use -echo screen as a -command-line option when running LAMMPS to see the offending line. - - */ diff --git a/src/GRANULAR/pair_granular_multi.cpp b/src/GRANULAR/pair_granular_multi.cpp deleted file mode 100644 index 07a6cab3dd..0000000000 --- a/src/GRANULAR/pair_granular_multi.cpp +++ /dev/null @@ -1,1624 +0,0 @@ -/* ---------------------------------------------------------------------- -http://lammps.sandia.gov, Sandia National Laboratories -Steve Plimpton, sjplimp@sandia.gov - -Copyright (2003) Sandia Corporation. Under the terms of Contract -DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains -certain rights in this software. This software is distributed under -the GNU General Public License. - -See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- -Contributing authors: -Dan Bolintineanu (SNL), Ishan Srivastava (SNL), Jeremy Lechman(SNL) -Leo Silbert (SNL), Gary Grest (SNL) ------------------------------------------------------------------------ */ - -#include -#include -#include -#include -#include "pair_granular_multi.h" -#include "atom.h" -#include "atom_vec.h" -#include "domain.h" -#include "force.h" -#include "update.h" -#include "modify.h" -#include "fix.h" -#include "fix_neigh_history.h" -#include "comm.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "neigh_request.h" -#include "memory.h" -#include "error.h" -#include "math_const.h" - -using namespace LAMMPS_NS; -using namespace MathConst; - -#define PI27SQ 266.47931882941264802866 // 27*PI**2 -#define THREEROOT3 5.19615242270663202362 // 3*sqrt(3) -#define SIXROOT6 14.69693845669906728801 // 6*sqrt(6) -#define INVROOT6 0.40824829046386307274 // 1/sqrt(6) -#define FOURTHIRDS 1.333333333333333 // 4/3 -#define TWOPI 6.28318530717959 // 2*PI - -#define EPSILON 1e-10 - -enum {HOOKE, HERTZ, HERTZ_MATERIAL, DMT, JKR}; -enum {VELOCITY, VISCOELASTIC, TSUJI}; -enum {TANGENTIAL_NOHISTORY, TANGENTIAL_MINDLIN}; -enum {TWIST_NONE, TWIST_NOHISTORY, TWIST_SDS, TWIST_MARSHALL}; -enum {ROLL_NONE, ROLL_NOHISTORY, ROLL_SDS}; - -/* ---------------------------------------------------------------------- */ - -PairGranularMulti::PairGranularMulti(LAMMPS *lmp) : Pair(lmp) -{ - single_enable = 1; - no_virial_fdotr_compute = 1; - fix_history = NULL; - - single_extra = 9; - svector = new double[single_extra]; - - neighprev = 0; - - nmax = 0; - mass_rigid = NULL; - - onerad_dynamic = NULL; - onerad_frozen = NULL; - maxrad_dynamic = NULL; - maxrad_frozen = NULL; - - dt = update->dt; - - // set comm size needed by this Pair if used with fix rigid - - comm_forward = 1; - - use_history = 0; - beyond_contact = 0; - nondefault_history_transfer = 0; - tangential_history_index = 0; - roll_history_index = twist_history_index = 0; - -} - -/* ---------------------------------------------------------------------- */ -PairGranularMulti::~PairGranularMulti() -{ - delete [] svector; - if (fix_history) modify->delete_fix("NEIGH_HISTORY"); - - if (allocated) { - memory->destroy(setflag); - memory->destroy(cutsq); - memory->destroy(cut); - - memory->destroy(normal_coeffs); - memory->destroy(tangential_coeffs); - memory->destroy(roll_coeffs); - memory->destroy(twist_coeffs); - - memory->destroy(normal); - memory->destroy(damping); - memory->destroy(tangential); - memory->destroy(roll); - memory->destroy(twist); - - delete [] onerad_dynamic; - delete [] onerad_frozen; - delete [] maxrad_dynamic; - delete [] maxrad_frozen; - } - memory->destroy(mass_rigid); -} - -void PairGranularMulti::compute(int eflag, int vflag) -{ - int i,j,ii,jj,inum,jnum,itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,nx,ny,nz; - double radi,radj,radsum,rsq,r,rinv,rsqinv; - double Reff, delta, dR, dR2; - - double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; - double wr1,wr2,wr3; - double vtr1,vtr2,vtr3,vrel; - - double knfac, damp_normal; - double k_tangential, damp_tangential; - double Fne, Ft, Fdamp, Fntot, Fcrit, Fscrit, Frcrit; - double fs, fs1, fs2, fs3; - - //For JKR - double R2, coh, F_pulloff, delta_pulloff, dist_pulloff, a, a2, E; - double t0, t1, t2, t3, t4, t5, t6; - double sqrt1, sqrt2, sqrt3, sqrt4; - - double mi,mj,meff,damp,ccel,tor1,tor2,tor3; - double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3,vrlmag,vrlmaginv; - - //Rolling - double k_roll, damp_roll; - double roll1, roll2, roll3, torroll1, torroll2, torroll3; - double rollmag, rolldotn, scalefac; - double fr, fr1, fr2, fr3; - - //Twisting - double k_twist, damp_twist, mu_twist; - double signtwist, magtwist, magtortwist, Mtcrit; - double tortwist1, tortwist2, tortwist3; - - double shrmag,rsht; - int *ilist,*jlist,*numneigh,**firstneigh; - int *touch,**firsttouch; - double *history,*allhistory,**firsthistory; - - bool touchflag; - - if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = vflag_fdotr = 0; - - int historyupdate = 1; - if (update->setupflag) historyupdate = 0; - - // update rigid body info for owned & ghost atoms if using FixRigid masses - // body[i] = which body atom I is in, -1 if none - // mass_body = mass of each rigid body - - if (fix_rigid && neighbor->ago == 0){ - int tmp; - int *body = (int *) fix_rigid->extract("body",tmp); - double *mass_body = (double *) fix_rigid->extract("masstotal",tmp); - if (atom->nmax > nmax) { - memory->destroy(mass_rigid); - nmax = atom->nmax; - memory->create(mass_rigid,nmax,"pair:mass_rigid"); - } - int nlocal = atom->nlocal; - for (i = 0; i < nlocal; i++) - if (body[i] >= 0) mass_rigid[i] = mass_body[body[i]]; - else mass_rigid[i] = 0.0; - comm->forward_comm_pair(this); - } - - double **x = atom->x; - double **v = atom->v; - double **f = atom->f; - int *type = atom->type; - double **omega = atom->omega; - double **torque = atom->torque; - double *radius = atom->radius; - double *rmass = atom->rmass; - int *mask = atom->mask; - int nlocal = atom->nlocal; - int newton_pair = force->newton_pair; - - inum = list->inum; - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->firstneigh; - firsttouch = fix_history->firstflag; - firsthistory = fix_history->firstvalue; - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - itype = type[i]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - itype = type[i]; - radi = radius[i]; - touch = firsttouch[i]; - allhistory = firsthistory[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for (jj = 0; jj < jnum; jj++){ - j = jlist[jj]; - j &= NEIGHMASK; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - jtype = type[j]; - rsq = delx*delx + dely*dely + delz*delz; - radj = radius[j]; - radsum = radi + radj; - - E = normal_coeffs[itype][jtype][0]; - Reff = radi*radj/(radi+radj); - touchflag = false; - - if (normal[itype][jtype] == JKR){ - if (touch[jj]){ - R2 = Reff*Reff; - coh = normal_coeffs[itype][jtype][3]; - a = cbrt(9.0*M_PI*coh*R2/(4*E)); - delta_pulloff = a*a/Reff - 2*sqrt(M_PI*coh*a/E); - dist_pulloff = radsum-delta_pulloff; - touchflag = (rsq < dist_pulloff*dist_pulloff); - } - else{ - touchflag = (rsq < radsum*radsum); - } - } - else{ - touchflag = (rsq < radsum*radsum); - } - - if (!touchflag){ - // unset non-touching neighbors - touch[jj] = 0; - history = &allhistory[size_history*jj]; - for (int k = 0; k < size_history; k++) history[k] = 0.0; - } - else{ - r = sqrt(rsq); - rinv = 1.0/r; - - nx = delx*rinv; - ny = dely*rinv; - nz = delz*rinv; - - // relative translational velocity - - vr1 = v[i][0] - v[j][0]; - vr2 = v[i][1] - v[j][1]; - vr3 = v[i][2] - v[j][2]; - - // normal component - - vnnr = vr1*nx + vr2*ny + vr3*nz; //v_R . n - vn1 = nx*vnnr; - vn2 = ny*vnnr; - vn3 = nz*vnnr; - - // meff = effective mass of pair of particles - // if I or J part of rigid body, use body mass - // if I or J is frozen, meff is other particle - - mi = rmass[i]; - mj = rmass[j]; - if (fix_rigid) { - if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; - if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; - } - - meff = mi*mj / (mi+mj); - if (mask[i] & freeze_group_bit) meff = mj; - if (mask[j] & freeze_group_bit) meff = mi; - - delta = radsum - r; - dR = delta*Reff; - if (normal[itype][jtype] == JKR){ - touch[jj] = 1; - R2=Reff*Reff; - coh = normal_coeffs[itype][jtype][3]; - dR2 = dR*dR; - t0 = coh*coh*R2*R2*E; - t1 = PI27SQ*t0; - t2 = 8*dR*dR2*E*E*E; - t3 = 4*dR2*E; - sqrt1 = MAX(0, t0*(t1+2*t2)); //In case of sqrt(0) < 0 due to precision issues - t4 = cbrt(t1+t2+THREEROOT3*M_PI*sqrt(sqrt1)); - t5 = t3/t4 + t4/E; - sqrt2 = MAX(0, 2*dR + t5); - t6 = sqrt(sqrt2); - sqrt3 = MAX(0, 4*dR - t5 + SIXROOT6*coh*M_PI*R2/(E*t6)); - a = INVROOT6*(t6 + sqrt(sqrt3)); - a2 = a*a; - knfac = FOURTHIRDS*E*a; - Fne = knfac*a2/Reff - TWOPI*a2*sqrt(4*coh*E/(M_PI*a)); - } - else{ - knfac = E; //Hooke - Fne = knfac*delta; - if (normal[itype][jtype] != HOOKE) - a = sqrt(dR); - Fne *= a; - if (normal[itype][jtype] == DMT) - Fne -= 4*MY_PI*normal_coeffs[itype][jtype][3]*Reff; - } - - //Consider restricting Hooke to only have 'velocity' as an option for damping? - if (damping[itype][jtype] == VELOCITY){ - damp_normal = 1; - } - else if (damping[itype][jtype] == VISCOELASTIC){ - if (normal[itype][jtype] == HOOKE) a = sqrt(dR); - damp_normal = a*meff; - } - else if (damping[itype][jtype] == TSUJI){ - damp_normal = sqrt(meff*knfac); - } - - Fdamp = -normal_coeffs[itype][jtype][1]*damp_normal*vnnr; - - Fntot = Fne + Fdamp; - - //**************************************** - //Tangential force, including history effects - //**************************************** - - // tangential component - vt1 = vr1 - vn1; - vt2 = vr2 - vn2; - vt3 = vr3 - vn3; - - // relative rotational velocity - wr1 = (radi*omega[i][0] + radj*omega[j][0]); - wr2 = (radi*omega[i][1] + radj*omega[j][1]); - wr3 = (radi*omega[i][2] + radj*omega[j][2]); - - // relative tangential velocities - vtr1 = vt1 - (nz*wr2-ny*wr3); - vtr2 = vt2 - (nx*wr3-nz*wr1); - vtr3 = vt3 - (ny*wr1-nx*wr2); - vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; - vrel = sqrt(vrel); - - // If any history is needed: - if (use_history){ - touch[jj] = 1; - history = &allhistory[size_history*jj]; - } - - if (normal[itype][jtype] == JKR){ - F_pulloff = 3*M_PI*coh*Reff; - Fcrit = fabs(Fne + 2*F_pulloff); - } - else{ - Fcrit = fabs(Fne); - } - - //------------------------------ - //Tangential forces - //------------------------------ - k_tangential = tangential_coeffs[itype][jtype][0]; - damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal; - - if (tangential_history){ - shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + - history[2]*history[2]); - - // Rotate and update displacements. - // See e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 - if (historyupdate) { - rsht = history[0]*nx + history[1]*ny + history[2]*nz; - if (fabs(rsht) < EPSILON) rsht = 0; - if (rsht > 0){ - scalefac = shrmag/(shrmag - rsht); //if rhst == shrmag, contacting pair has rotated 90 deg. in one step, in which case you deserve a crash! - history[0] -= rsht*nx; - history[1] -= rsht*ny; - history[2] -= rsht*nz; - //Also rescale to preserve magnitude - history[0] *= scalefac; - history[1] *= scalefac; - history[2] *= scalefac; - } - //Update history - history[0] += vtr1*dt; - history[1] += vtr2*dt; - history[2] += vtr3*dt; - } - - // tangential forces = history + tangential velocity damping - fs1 = -k_tangential*history[0] - damp_tangential*vtr1; - fs2 = -k_tangential*history[1] - damp_tangential*vtr2; - fs3 = -k_tangential*history[2] - damp_tangential*vtr3; - - // rescale frictional displacements and forces if needed - Fscrit = tangential_coeffs[itype][jtype][2] * Fcrit; - fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); - if (fs > Fscrit) { - if (shrmag != 0.0) { - history[0] = -1.0/k_tangential*(Fscrit*fs1/fs + damp_tangential*vtr1); - history[1] = -1.0/k_tangential*(Fscrit*fs2/fs + damp_tangential*vtr2); - history[2] = -1.0/k_tangential*(Fscrit*fs3/fs + damp_tangential*vtr3); - fs1 *= Fscrit/fs; - fs2 *= Fscrit/fs; - fs3 *= Fscrit/fs; - } else fs1 = fs2 = fs3 = 0.0; - } - } - else{ //Classic pair gran/hooke (no history) - fs = meff*damp_tangential*vrel; - if (vrel != 0.0) Ft = MIN(Fne,fs) / vrel; - else Ft = 0.0; - fs1 = -Ft*vtr1; - fs2 = -Ft*vtr2; - fs3 = -Ft*vtr3; - } - - //**************************************** - // Rolling resistance - //**************************************** - - if (roll[itype][jtype] != ROLL_NONE){ - relrot1 = omega[i][0] - omega[j][0]; - relrot2 = omega[i][1] - omega[j][1]; - relrot3 = omega[i][2] - omega[j][2]; - - // rolling velocity, see eq. 31 of Wang et al, Particuology v 23, p 49 (2015) - // This is different from the Marshall papers, which use the Bagi/Kuhn formulation - // for rolling velocity (see Wang et al for why the latter is wrong) - vrl1 = Reff*(relrot2*nz - relrot3*ny); //- 0.5*((radj-radi)/radsum)*vtr1; - vrl2 = Reff*(relrot3*nx - relrot1*nz); //- 0.5*((radj-radi)/radsum)*vtr2; - vrl3 = Reff*(relrot1*ny - relrot2*nx); //- 0.5*((radj-radi)/radsum)*vtr3; - vrlmag = sqrt(vrl1*vrl1+vrl2*vrl2+vrl3*vrl3); - if (vrlmag != 0.0) vrlmaginv = 1.0/vrlmag; - else vrlmaginv = 0.0; - - if (roll_history){ - int rhist0 = roll_history_index; - int rhist1 = rhist0 + 1; - int rhist2 = rhist1 + 1; - - // Rolling displacement - rollmag = sqrt(history[rhist0]*history[rhist0] + - history[rhist1]*history[rhist1] + - history[rhist2]*history[rhist2]); - - rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz; - - if (historyupdate){ - if (fabs(rolldotn) < EPSILON) rolldotn = 0; - if (rolldotn > 0){ //Rotate into tangential plane - scalefac = rollmag/(rollmag - rolldotn); - history[rhist0] -= rolldotn*nx; - history[rhist1] -= rolldotn*ny; - history[rhist2] -= rolldotn*nz; - //Also rescale to preserve magnitude - history[rhist0] *= scalefac; - history[rhist1] *= scalefac; - history[rhist2] *= scalefac; - } - history[rhist0] += vrl1*dt; - history[rhist1] += vrl2*dt; - history[rhist2] += vrl3*dt; - } - - - k_roll = roll_coeffs[itype][jtype][0]; - damp_roll = roll_coeffs[itype][jtype][1]; - fr1 = -k_roll*history[rhist0] - damp_roll*vrl1; - fr2 = -k_roll*history[rhist1] - damp_roll*vrl2; - fr3 = -k_roll*history[rhist2] - damp_roll*vrl3; - - // rescale frictional displacements and forces if needed - Frcrit = roll_coeffs[itype][jtype][2] * Fcrit; - - fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3); - if (fr > Frcrit) { - if (rollmag != 0.0) { - history[rhist0] = -1.0/k_roll*(Frcrit*fr1/fr + damp_roll*vrl1); - history[rhist1] = -1.0/k_roll*(Frcrit*fr2/fr + damp_roll*vrl2); - history[rhist2] = -1.0/k_roll*(Frcrit*fr3/fr + damp_roll*vrl3); - fr1 *= Frcrit/fr; - fr2 *= Frcrit/fr; - fr3 *= Frcrit/fr; - } else fr1 = fr2 = fr3 = 0.0; - } - } - else{ // - fr = meff*roll_coeffs[itype][jtype][1]*vrlmag; - if (vrlmag != 0.0) fr = MIN(Fne, fr) / vrlmag; - else fr = 0.0; - fr1 = -fr*vrl1; - fr2 = -fr*vrl2; - fr3 = -fr*vrl3; - } - } - - //**************************************** - // Twisting torque, including history effects - //**************************************** - if (twist[itype][jtype] != TWIST_NONE){ - magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of Marshall) - if (twist[itype][jtype] == TWIST_MARSHALL){ - k_twist = 0.5*k_tangential*a*a;; //eq 32 - damp_twist = 0.5*damp_tangential*a*a; - mu_twist = TWOTHIRDS*a; - } - else{ - k_twist = twist_coeffs[itype][jtype][0]; - damp_twist = twist_coeffs[itype][jtype][1]; - mu_twist = twist_coeffs[itype][jtype][2]; - } - if (twist[itype][jtype] > 1){ - if (historyupdate){ - history[twist_history_index] += magtwist*dt; - } - magtortwist = -k_twist*history[twist_history_index] - damp_twist*magtwist;//M_t torque (eq 30) - signtwist = (magtwist > 0) - (magtwist < 0); - Mtcrit = TWOTHIRDS*a*Fscrit;//critical torque (eq 44) - if (fabs(magtortwist) > Mtcrit) { - history[twist_history_index] = 1.0/k_twist*(Mtcrit*signtwist - damp_twist*magtwist); - magtortwist = -Mtcrit * signtwist; //eq 34 - } - } - else{ - if (magtwist > 0) magtortwist = -damp_twist*magtwist; - else magtortwist = 0; - } - } - // Apply forces & torques - - fx = nx*Fntot + fs1; - fy = ny*Fntot + fs2; - fz = nz*Fntot + fs3; - - f[i][0] += fx; - f[i][1] += fy; - f[i][2] += fz; - - tor1 = ny*fs3 - nz*fs2; - tor2 = nz*fs1 - nx*fs3; - tor3 = nx*fs2 - ny*fs1; - - torque[i][0] -= radi*tor1; - torque[i][1] -= radi*tor2; - torque[i][2] -= radi*tor3; - - if (twist[itype][jtype] != TWIST_NONE){ - tortwist1 = magtortwist * nx; - tortwist2 = magtortwist * ny; - tortwist3 = magtortwist * nz; - - torque[i][0] += tortwist1; - torque[i][1] += tortwist2; - torque[i][2] += tortwist3; - } - - if (roll[itype][jtype] != ROLL_NONE){ - torroll1 = Reff*(ny*fr3 - nz*fr2); //n cross fr - torroll2 = Reff*(nz*fr1 - nx*fr3); - torroll3 = Reff*(nx*fr2 - ny*fr1); - - torque[i][0] += torroll1; - torque[i][1] += torroll2; - torque[i][2] += torroll3; - } - - if (force->newton_pair || j < nlocal) { - f[j][0] -= fx; - f[j][1] -= fy; - f[j][2] -= fz; - - torque[j][0] -= radj*tor1; - torque[j][1] -= radj*tor2; - torque[j][2] -= radj*tor3; - - if (twist[itype][jtype] != TWIST_NONE){ - torque[j][0] -= tortwist1; - torque[j][1] -= tortwist2; - torque[j][2] -= tortwist3; - } - if (roll[itype][jtype] != ROLL_NONE){ - torque[j][0] -= torroll1; - torque[j][1] -= torroll2; - torque[j][2] -= torroll3; - } - } - if (evflag) ev_tally_xyz(i,j,nlocal,0, - 0.0,0.0,fx,fy,fz,delx,dely,delz); - } - } - } -} - - -/* ---------------------------------------------------------------------- -allocate all arrays -------------------------------------------------------------------------- */ - -void PairGranularMulti::allocate() -{ - allocated = 1; - int n = atom->ntypes; - - memory->create(setflag,n+1,n+1,"pair:setflag"); - for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) - setflag[i][j] = 0; - - memory->create(cutsq,n+1,n+1,"pair:cutsq"); - memory->create(cut,n+1,n+1,"pair:cut"); - memory->create(normal_coeffs,n+1,n+1,4,"pair:normal_coeffs"); - memory->create(tangential_coeffs,n+1,n+1,3,"pair:tangential_coeffs"); - memory->create(roll_coeffs,n+1,n+1,3,"pair:roll_coeffs"); - memory->create(twist_coeffs,n+1,n+1,3,"pair:twist_coeffs"); - - memory->create(normal,n+1,n+1,"pair:normal"); - memory->create(damping,n+1,n+1,"pair:damping"); - memory->create(tangential,n+1,n+1,"pair:tangential"); - memory->create(roll,n+1,n+1,"pair:roll"); - memory->create(twist,n+1,n+1,"pair:twist"); - - onerad_dynamic = new double[n+1]; - onerad_frozen = new double[n+1]; - maxrad_dynamic = new double[n+1]; - maxrad_frozen = new double[n+1]; -} - -/* ---------------------------------------------------------------------- - global settings -------------------------------------------------------------------------- */ - -void PairGranularMulti::settings(int narg, char **arg) -{ - if (narg == 1){ - cutoff_global = force->numeric(FLERR,arg[0]); - } - else{ - cutoff_global = -1; //Will be set based on particle sizes, model choice - } - - tangential_history = 0; - roll_history = twist_history = 0; -} - -/* ---------------------------------------------------------------------- - set coeffs for one or more type pairs -------------------------------------------------------------------------- */ - -void PairGranularMulti::coeff(int narg, char **arg) -{ - int normal_local, damping_local, tangential_local, roll_local, twist_local; - - double normal_coeffs_local[4]; - double tangential_coeffs_local[4]; - double roll_coeffs_local[4]; - double twist_coeffs_local[4]; - - if (narg < 2) - error->all(FLERR,"Incorrect args for pair coefficients"); - - if (!allocated) allocate(); - - int ilo,ihi,jlo,jhi; - force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi); - force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi); - - //Defaults - normal_local = tangential_local = -1; - roll_local = twist_local = 0; - damping_local = VISCOELASTIC; - - int iarg = 2; - while (iarg < narg){ - if (strcmp(arg[iarg], "hooke") == 0){ - if (iarg + 2 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hooke option"); - normal_local = HOOKE; - normal_coeffs_local[0] = force->numeric(FLERR,arg[iarg+1]); //kn - normal_coeffs_local[1] = force->numeric(FLERR,arg[iarg+2]); //damping - iarg += 3; - } - else if (strcmp(arg[iarg], "hertz") == 0){ - int num_coeffs = 2; - if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hertz option"); - normal_local = HERTZ; - normal_coeffs_local[0] = force->numeric(FLERR,arg[iarg+1]); //kn - normal_coeffs_local[1] = force->numeric(FLERR,arg[iarg+2]); //damping - iarg += num_coeffs+1; - } - else if (strcmp(arg[iarg], "hertz/material") == 0){ - int num_coeffs = 3; - if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hertz option"); - normal_local = HERTZ; - normal_coeffs_local[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //E - normal_coeffs_local[1] = force->numeric(FLERR,arg[iarg+2]); //damping - normal_coeffs_local[2] = force->numeric(FLERR,arg[iarg+3]); //G - iarg += num_coeffs+1; - } - else if (strcmp(arg[iarg], "dmt") == 0){ - if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hertz option"); - normal_local = DMT; - normal_coeffs_local[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //E - normal_coeffs_local[1] = force->numeric(FLERR,arg[iarg+2]); //damping - normal_coeffs_local[2] = force->numeric(FLERR,arg[iarg+3]); //G - normal_coeffs_local[3] = force->numeric(FLERR,arg[iarg+3]); //cohesion - iarg += 5; - } - else if (strcmp(arg[iarg], "jkr") == 0){ - if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for JKR option"); - beyond_contact = 1; - normal_local = JKR; - normal_coeffs_local[0] = force->numeric(FLERR,arg[iarg+1]); //E - normal_coeffs_local[1] = force->numeric(FLERR,arg[iarg+2]); //damping - normal_coeffs_local[2] = force->numeric(FLERR,arg[iarg+3]); //G - normal_coeffs_local[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion - iarg += 5; - } - else if (strcmp(arg[iarg], "damp") == 0){ - if (iarg+1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters provided for damping model"); - if (strcmp(arg[iarg+1], "velocity") == 0){ - damping_local = VELOCITY; - iarg += 1; - } - else if (strcmp(arg[iarg+1], "viscoelastic") == 0){ - damping_local = VISCOELASTIC; - iarg += 1; - } - else if (strcmp(arg[iarg], "tsuji") == 0){ - damping_local = TSUJI; - iarg += 1; - } - } - else if (strcmp(arg[iarg], "tangential") == 0){ - if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model"); - if (strcmp(arg[iarg+1], "nohistory") == 0){ - tangential_local = TANGENTIAL_NOHISTORY; - } - else if (strcmp(arg[iarg+1], "mindlin") == 0){ - tangential_local = TANGENTIAL_MINDLIN; - tangential_history = 1; - } - else{ - error->all(FLERR, "Illegal pair_coeff command, tangential model not recognized"); - } - tangential_coeffs_local[0] = force->numeric(FLERR,arg[iarg+2]); //kt - tangential_coeffs_local[1] = force->numeric(FLERR,arg[iarg+3]); //gammat - tangential_coeffs_local[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. - iarg += 5; - } - else if (strcmp(arg[iarg], "rolling") == 0){ - if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); - if (strcmp(arg[iarg+1], "none") == 0){ - roll_local = ROLL_NONE; - iarg += 2; - } - else{ - if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for rolling model"); - if (strcmp(arg[iarg+1], "nohistory") == 0){ - roll_local = ROLL_NOHISTORY; - } - else if (strcmp(arg[iarg+1], "sds") == 0){ - roll_local = ROLL_SDS; - roll_history = 1; - } - else{ - error->all(FLERR, "Illegal pair_coeff command, rolling friction model not recognized"); - } - roll_coeffs_local[0] = force->numeric(FLERR,arg[iarg+2]); //kR - roll_coeffs_local[1] = force->numeric(FLERR,arg[iarg+3]); //gammaR - roll_coeffs_local[2] = force->numeric(FLERR,arg[iarg+4]); //rolling friction coeff. - iarg += 5; - } - } - else if (strcmp(arg[iarg], "twisting") == 0){ - if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); - if (strcmp(arg[iarg+1], "none") == 0){ - twist_local = TWIST_NONE; - iarg += 2; - } - else if (strcmp(arg[iarg+1], "marshall") == 0){ - twist_local = TWIST_MARSHALL; - twist_history = 1; - iarg += 2; - } - else{ - if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for twist model"); - if (strcmp(arg[iarg+1], "nohistory") == 0){ - twist_local = TWIST_NOHISTORY; - } - else if (strcmp(arg[iarg+1], "sds") == 0){ - twist_local = TWIST_SDS; - twist_history = 1; - } - else{ - error->all(FLERR, "Illegal pair_coeff command, twisting friction model not recognized"); - } - twist_coeffs_local[0] = force->numeric(FLERR,arg[iarg+2]); //kt - twist_coeffs_local[1] = force->numeric(FLERR,arg[iarg+3]); //gammat - twist_coeffs_local[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. - iarg += 5; - } - } - else error->all(FLERR, "Illegal pair coeff command"); - } - - //It is an error not to specify normal or tangential model - if ((normal_local < 0) || (tangential_local < 0)) error->all(FLERR, "Illegal pair coeff command, must specify normal contact model"); - - int count = 0; - double damp; - if (damping_local == TSUJI){ - double cor; - cor = normal_coeffs_local[1]; - damp = 1.2728-4.2783*cor+11.087*pow(cor,2)-22.348*pow(cor,3)+ - 27.467*pow(cor,4)-18.022*pow(cor,5)+ - 4.8218*pow(cor,6); - } - else damp = normal_coeffs_local[1]; - - for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { - normal[i][j] = normal[j][i] = normal_local; - normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = normal_coeffs_local[0]; - normal_coeffs[i][j][1] = normal_coeffs[j][i][1] = damp; - if (normal_local != HERTZ && normal_local != HOOKE) normal_coeffs[i][j][2] = normal_coeffs_local[2]; - if ((normal_local == JKR) || (normal_local == DMT)) - normal_coeffs[i][j][3] = normal_coeffs[j][i][3] = normal_coeffs_local[3]; - - damping[i][j] = damping[j][i] = damping_local; - - tangential[i][j] = tangential[j][i] = tangential_local; - for (int k = 0; k < 3; k++) - tangential_coeffs[i][j][k] = tangential_coeffs[j][i][k] = tangential_coeffs_local[k]; - - roll[i][j] = roll[j][i] = roll_local; - if (roll_local != ROLL_NONE) - for (int k = 0; k < 3; k++) - roll_coeffs[i][j][k] = roll_coeffs[j][i][k] = roll_coeffs_local[k]; - - twist[i][j] = twist[j][i] = twist_local; - if (twist_local != TWIST_NONE && twist_local != TWIST_MARSHALL) - for (int k = 0; k < 3; k++) - twist_coeffs[i][j][k] = twist_coeffs[j][i][k] = twist_coeffs_local[k]; - - setflag[i][j] = 1; - count++; - } - } - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); -} - -/* ---------------------------------------------------------------------- - init specific to this pair style -------------------------------------------------------------------------- */ - -void PairGranularMulti::init_style() -{ - int i; - - // error and warning checks - - if (!atom->radius_flag || !atom->rmass_flag) - error->all(FLERR,"Pair granular requires atom attributes radius, rmass"); - if (comm->ghost_velocity == 0) - error->all(FLERR,"Pair granular requires ghost atoms store velocity"); - - // Determine whether we need a granular neigh list, how large it needs to be - use_history = tangential_history || roll_history || twist_history; - - //For JKR, will need fix/neigh/history to keep track of touch arrays - for (int i = 1; i <= atom->ntypes; i++) - for (int j = 1; j <= atom->ntypes; j++) - if (normal[i][j] == JKR) use_history = 1; - - size_history = 3*tangential_history + 3*roll_history + twist_history; - - //Determine location of tangential/roll/twist histories in array - if (roll_history){ - if (tangential_history) roll_history_index = 3; - else roll_history_index = 0; - } - if (twist_history){ - if (tangential_history){ - if (roll_history) twist_history_index = 6; - else twist_history_index = 3; - } - else{ - if (roll_history) twist_history_index = 3; - else twist_history_index = 0; - } - } - - int irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->size = 1; - if (use_history) neighbor->requests[irequest]->history = 1; - - dt = update->dt; - - // if history is stored: - // if first init, create Fix needed for storing history - - if (use_history && fix_history == NULL) { - char dnumstr[16]; - sprintf(dnumstr,"%d",size_history); - char **fixarg = new char*[4]; - fixarg[0] = (char *) "NEIGH_HISTORY"; - fixarg[1] = (char *) "all"; - fixarg[2] = (char *) "NEIGH_HISTORY"; - fixarg[3] = dnumstr; - modify->add_fix(4,fixarg,1); - delete [] fixarg; - fix_history = (FixNeighHistory *) modify->fix[modify->nfix-1]; - fix_history->pair = this; - } - - // check for FixFreeze and set freeze_group_bit - - for (i = 0; i < modify->nfix; i++) - if (strcmp(modify->fix[i]->style,"freeze") == 0) break; - if (i < modify->nfix) freeze_group_bit = modify->fix[i]->groupbit; - else freeze_group_bit = 0; - - // check for FixRigid so can extract rigid body masses - - fix_rigid = NULL; - for (i = 0; i < modify->nfix; i++) - if (modify->fix[i]->rigid_flag) break; - if (i < modify->nfix) fix_rigid = modify->fix[i]; - - // check for FixPour and FixDeposit so can extract particle radii - - int ipour; - for (ipour = 0; ipour < modify->nfix; ipour++) - if (strcmp(modify->fix[ipour]->style,"pour") == 0) break; - if (ipour == modify->nfix) ipour = -1; - - int idep; - for (idep = 0; idep < modify->nfix; idep++) - if (strcmp(modify->fix[idep]->style,"deposit") == 0) break; - if (idep == modify->nfix) idep = -1; - - // set maxrad_dynamic and maxrad_frozen for each type - // include future FixPour and FixDeposit particles as dynamic - - int itype; - for (i = 1; i <= atom->ntypes; i++) { - onerad_dynamic[i] = onerad_frozen[i] = 0.0; - if (ipour >= 0) { - itype = i; - double radmax = *((double *) modify->fix[ipour]->extract("radius",itype)); - if (normal[itype][itype] == JKR) radmax = radmax - 0.5*pulloff_distance(radmax, itype); - onerad_dynamic[i] = radmax; - } - if (idep >= 0) { - itype = i; - double radmax = *((double *) modify->fix[idep]->extract("radius",itype)); - if (normal[itype][itype] == JKR) radmax = radmax - 0.5*pulloff_distance(radmax, itype); - onerad_dynamic[i] = radmax; - } - } - - double *radius = atom->radius; - int *mask = atom->mask; - int *type = atom->type; - int nlocal = atom->nlocal; - - for (i = 0; i < nlocal; i++){ - double radius_cut = radius[i]; - if (normal[type[i]][type[i]] == JKR){ - radius_cut = radius[i] - 0.5*pulloff_distance(radius[i], type[i]); - } - if (mask[i] & freeze_group_bit){ - onerad_frozen[type[i]] = MAX(onerad_frozen[type[i]],radius_cut); - } - else{ - onerad_dynamic[type[i]] = MAX(onerad_dynamic[type[i]],radius_cut); - } - } - - MPI_Allreduce(&onerad_dynamic[1],&maxrad_dynamic[1],atom->ntypes, - MPI_DOUBLE,MPI_MAX,world); - MPI_Allreduce(&onerad_frozen[1],&maxrad_frozen[1],atom->ntypes, - MPI_DOUBLE,MPI_MAX,world); - - // set fix which stores history info - - if (size_history > 0){ - int ifix = modify->find_fix("NEIGH_HISTORY"); - if (ifix < 0) error->all(FLERR,"Could not find pair fix neigh history ID"); - fix_history = (FixNeighHistory *) modify->fix[ifix]; - } -} - -/* ---------------------------------------------------------------------- - init for one type pair i,j and corresponding j,i -------------------------------------------------------------------------- */ - -double PairGranularMulti::init_one(int i, int j) -{ - double cutoff; - if (setflag[i][j] == 0) { - if ((normal[i][i] != normal[j][j]) || - (damping[i][i] != damping[j][j]) || - (tangential[i][i] != tangential[j][j]) || - (roll[i][i] != roll[j][j]) || - (twist[i][i] != twist[j][j])){ - - char str[512]; - sprintf(str,"Granular pair style functional forms are different, cannot mix coefficients for types %d and %d. \nThis combination must be set explicitly via pair_coeff command.",i,j); - error->one(FLERR,str); - } - - if (normal[i][j] != HOOKE && normal[i][j] != HERTZ){ - normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = mix_stiffnessE(normal_coeffs[i][i][0], normal_coeffs[j][j][0], - normal_coeffs[i][i][2], normal_coeffs[j][j][2]); - normal_coeffs[i][j][2] = normal_coeffs[j][i][2] = mix_stiffnessG(normal_coeffs[i][i][0], normal_coeffs[j][j][0], - normal_coeffs[i][i][2], normal_coeffs[j][j][2]); - } - else{ - normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = mix_geom(normal_coeffs[i][i][0], normal_coeffs[j][j][0]); - } - - normal_coeffs[i][j][1] = normal_coeffs[j][i][1] = mix_geom(normal_coeffs[i][i][1], normal_coeffs[j][j][1]); - if ((normal[i][i] == JKR) || (normal[i][i] == DMT)) - normal_coeffs[i][j][3] = normal_coeffs[j][i][3] = mix_geom(normal_coeffs[i][i][3], normal_coeffs[j][j][3]); - - for (int k = 0; k < 3; k++) - tangential_coeffs[i][j][k] = normal_coeffs[j][i][k] = mix_geom(tangential_coeffs[i][i][k], tangential_coeffs[j][j][k]); - - - if (roll[i][i] != ROLL_NONE){ - for (int k = 0; k < 3; k++) - roll_coeffs[i][j][k] = roll_coeffs[j][i][k] = mix_geom(roll_coeffs[i][i][k], roll_coeffs[j][j][k]); - } - - if (twist[i][i] != TWIST_NONE && twist[i][i] != TWIST_MARSHALL){ - for (int k = 0; k < 3; k++) - twist_coeffs[i][j][k] = twist_coeffs[j][i][k] = mix_geom(twist_coeffs[i][i][k], twist_coeffs[j][j][k]); - } - } - - // It is possible that cut[i][j] at this point is still 0.0. This can happen when - // there is a future fix_pour after the current run. A cut[i][j] = 0.0 creates - // problems because neighbor.cpp uses min(cut[i][j]) to decide on the bin size - // To avoid this issue, for cases involving cut[i][j] = 0.0 (possible only - // if there is no current information about radius/cutoff of type i and j). - // we assign cutoff = max(cut[i][j]) for i,j such that cut[i][j] > 0.0. - - if (cutoff_global < 0){ - if (((maxrad_dynamic[i] > 0.0) && (maxrad_dynamic[j] > 0.0)) || - ((maxrad_dynamic[i] > 0.0) && (maxrad_frozen[j] > 0.0)) || - ((maxrad_frozen[i] > 0.0) && (maxrad_dynamic[j] > 0.0))) { // radius info about both i and j exist - cutoff = maxrad_dynamic[i]+maxrad_dynamic[j]; - cutoff = MAX(cutoff,maxrad_frozen[i]+maxrad_dynamic[j]); - cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]); - } - else { // radius info about either i or j does not exist (i.e. not present and not about to get poured; set to largest value to not interfere with neighbor list) - double cutmax = 0.0; - for (int k = 1; k <= atom->ntypes; k++) { - cutmax = MAX(cutmax,2.0*maxrad_dynamic[k]); - cutmax = MAX(cutmax,2.0*maxrad_frozen[k]); - } - cutoff = cutmax; - } - } - else{ - cutoff = cutoff_global; - } - return cutoff; -} - - -/* ---------------------------------------------------------------------- - proc 0 writes to restart file - ------------------------------------------------------------------------- */ - -void PairGranularMulti::write_restart(FILE *fp) -{ - int i,j; - for (i = 1; i <= atom->ntypes; i++) { - for (j = i; j <= atom->ntypes; j++) { - fwrite(&setflag[i][j],sizeof(int),1,fp); - if (setflag[i][j]) { - fwrite(&normal[i][j],sizeof(int),1,fp); - fwrite(&damping[i][j],sizeof(int),1,fp); - fwrite(&tangential[i][j],sizeof(int),1,fp); - fwrite(&roll[i][j],sizeof(int),1,fp); - fwrite(&twist[i][j],sizeof(int),1,fp); - fwrite(&normal_coeffs[i][j],sizeof(double),4,fp); - fwrite(&tangential_coeffs[i][j],sizeof(double),3,fp); - fwrite(&roll_coeffs[i][j],sizeof(double),3,fp); - fwrite(&twist_coeffs[i][j],sizeof(double),3,fp); - fwrite(&cut[i][j],sizeof(double),1,fp); - } - } - } -} - -/* ---------------------------------------------------------------------- - proc 0 reads from restart file, bcasts - ------------------------------------------------------------------------- */ - -void PairGranularMulti::read_restart(FILE *fp) -{ - allocate(); - int i,j; - int me = comm->me; - for (i = 1; i <= atom->ntypes; i++) { - for (j = i; j <= atom->ntypes; j++) { - if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); - MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); - if (setflag[i][j]) { - if (me == 0) { - fread(&normal[i][j],sizeof(int),1,fp); - fread(&damping[i][j],sizeof(int),1,fp); - fread(&tangential[i][j],sizeof(int),1,fp); - fread(&roll[i][j],sizeof(int),1,fp); - fread(&twist[i][j],sizeof(int),1,fp); - fread(&normal_coeffs[i][j],sizeof(double),4,fp); - fread(&tangential_coeffs[i][j],sizeof(double),3,fp); - fread(&roll_coeffs[i][j],sizeof(double),3,fp); - fread(&twist_coeffs[i][j],sizeof(double),3,fp); - fread(&cut[i][j],sizeof(double),1,fp); - } - MPI_Bcast(&normal[i][j],1,MPI_INT,0,world); - MPI_Bcast(&damping[i][j],1,MPI_INT,0,world); - MPI_Bcast(&tangential[i][j],1,MPI_INT,0,world); - MPI_Bcast(&roll[i][j],1,MPI_INT,0,world); - MPI_Bcast(&twist[i][j],1,MPI_INT,0,world); - MPI_Bcast(&normal_coeffs[i][j],4,MPI_DOUBLE,0,world); - MPI_Bcast(&tangential_coeffs[i][j],3,MPI_DOUBLE,0,world); - MPI_Bcast(&roll_coeffs[i][j],3,MPI_DOUBLE,0,world); - MPI_Bcast(&twist_coeffs[i][j],3,MPI_DOUBLE,0,world); - MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); - } - } - } -} - - -/* ---------------------------------------------------------------------- */ - -void PairGranularMulti::reset_dt() -{ - dt = update->dt; -} - -/* ---------------------------------------------------------------------- */ - -double PairGranularMulti::single(int i, int j, int itype, int jtype, - double rsq, double factor_coul, double factor_lj, double &fforce) -{ - double radi,radj,radsum; - double r,rinv,rsqinv,delx,dely,delz, nx, ny, nz, Reff; - double dR, dR2; - double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3,wr1,wr2,wr3; - double vtr1,vtr2,vtr3,vrel; - double mi,mj,meff,damp,ccel,tor1,tor2,tor3; - double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3,vrlmag,vrlmaginv; - - double knfac, damp_normal; - double k_tangential, damp_tangential; - double Fne, Ft, Fdamp, Fntot, Fcrit, Fscrit, Frcrit; - double fs, fs1, fs2, fs3; - - //For JKR - double R2, coh, F_pulloff, delta_pulloff, dist_pulloff, a, a2, E; - double delta, t0, t1, t2, t3, t4, t5, t6; - double sqrt1, sqrt2, sqrt3, sqrt4; - - - //Rolling - double k_roll, damp_roll; - double roll1, roll2, roll3, torroll1, torroll2, torroll3; - double rollmag, rolldotn, scalefac; - double fr, fr1, fr2, fr3; - - //Twisting - double k_twist, damp_twist, mu_twist; - double signtwist, magtwist, magtortwist, Mtcrit; - double tortwist1, tortwist2, tortwist3; - - double shrmag,rsht; - int jnum; - int *ilist,*jlist,*numneigh,**firstneigh; - int *touch,**firsttouch; - double *history,*allhistory,**firsthistory; - - double *radius = atom->radius; - radi = radius[i]; - radj = radius[j]; - radsum = radi + radj; - Reff = radi*radj/(radi+radj); - - bool touchflag; - if (normal[itype][jtype] == JKR){ - R2 = Reff*Reff; - coh = normal_coeffs[itype][jtype][3]; - a = cbrt(9.0*M_PI*coh*R2/(4*E)); - delta_pulloff = a*a/Reff - 2*sqrt(M_PI*coh*a/E); - dist_pulloff = radsum+delta_pulloff; - touchflag = (rsq <= dist_pulloff*dist_pulloff); - } - else{ - touchflag = (rsq <= radsum*radsum); - } - - if (touchflag){ - fforce = 0.0; - for (int m = 0; m < single_extra; m++) svector[m] = 0.0; - return 0.0; - } - - double **x = atom->x; - delx = x[i][0] - x[j][0]; - dely = x[i][1] - x[j][1]; - delz = x[i][2] - x[j][2]; - r = sqrt(rsq); - rinv = 1.0/r; - - nx = delx*rinv; - ny = dely*rinv; - nz = delz*rinv; - - // relative translational velocity - - double **v = atom->v; - vr1 = v[i][0] - v[j][0]; - vr2 = v[i][1] - v[j][1]; - vr3 = v[i][2] - v[j][2]; - - // normal component - - vnnr = vr1*nx + vr2*ny + vr3*nz; - vn1 = nx*vnnr; - vn2 = ny*vnnr; - vn3 = nz*vnnr; - - double *rmass = atom->rmass; - int *mask = atom->mask; - mi = rmass[i]; - mj = rmass[j]; - if (fix_rigid) { - if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; - if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; - } - - meff = mi*mj / (mi+mj); - if (mask[i] & freeze_group_bit) meff = mj; - if (mask[j] & freeze_group_bit) meff = mi; - - delta = radsum - r; - dR = delta*Reff; - - // tangential component - - vt1 = vr1 - vn1; - vt2 = vr2 - vn2; - vt3 = vr3 - vn3; - - // relative rotational velocity - - double **omega = atom->omega; - wr1 = (radi*omega[i][0] + radj*omega[j][0]); - wr2 = (radi*omega[i][1] + radj*omega[j][1]); - wr3 = (radi*omega[i][2] + radj*omega[j][2]); - - // meff = effective mass of pair of particles - // if I or J part of rigid body, use body mass - // if I or J is frozen, meff is other particle - - int *type = atom->type; - - mi = rmass[i]; - mj = rmass[j]; - if (fix_rigid) { - // NOTE: ensure mass_rigid is current for owned+ghost atoms? - if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; - if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; - } - - meff = mi*mj / (mi+mj); - if (mask[i] & freeze_group_bit) meff = mj; - if (mask[j] & freeze_group_bit) meff = mi; - - delta = radsum - r; - dR = delta*Reff; - if (normal[itype][jtype] == JKR){ - dR2 = dR*dR; - t0 = coh*coh*R2*R2*E; - t1 = PI27SQ*t0; - t2 = 8*dR*dR2*E*E*E; - t3 = 4*dR2*E; - sqrt1 = MAX(0, t0*(t1+2*t2)); //In case of sqrt(0) < 0 due to precision issues - t4 = cbrt(t1+t2+THREEROOT3*M_PI*sqrt(sqrt1)); - t5 = t3/t4 + t4/E; - sqrt2 = MAX(0, 2*dR + t5); - t6 = sqrt(sqrt2); - sqrt3 = MAX(0, 4*dR - t5 + SIXROOT6*coh*M_PI*R2/(E*t6)); - a = INVROOT6*(t6 + sqrt(sqrt3)); - a2 = a*a; - knfac = FOURTHIRDS*E*a; - Fne = knfac*a2/Reff - TWOPI*a2*sqrt(4*coh*E/(M_PI*a)); - } - else{ - knfac = E; - Fne = knfac*delta; - if (normal[itype][jtype] != HOOKE) - a = sqrt(dR); - Fne *= a; - if (normal[itype][jtype] == DMT) - Fne -= 4*MY_PI*normal_coeffs[itype][jtype][3]*Reff; - } - - //Consider restricting Hooke to only have 'velocity' as an option for damping? - if (damping[itype][jtype] == VELOCITY){ - damp_normal = normal_coeffs[itype][jtype][1]; - } - else if (damping[itype][jtype] == VISCOELASTIC){ - if (normal[itype][jtype] == HOOKE) a = sqrt(dR); - damp_normal = normal_coeffs[itype][jtype][1]*a*meff; - } - else if (damping[itype][jtype] == TSUJI){ - damp_normal = normal_coeffs[itype][jtype][1]*sqrt(meff*knfac); - } - - Fdamp = -damp_normal*vnnr; - - Fntot = Fne + Fdamp; - - jnum = list->numneigh[i]; - jlist = list->firstneigh[i]; - - if (use_history){ - allhistory = fix_history->firstvalue[i]; - for (int jj = 0; jj < jnum; jj++) { - neighprev++; - if (neighprev >= jnum) neighprev = 0; - if (jlist[neighprev] == j) break; - } - history = &allhistory[size_history*neighprev]; - } - - //**************************************** - //Tangential force, including history effects - //**************************************** - - // tangential component - vt1 = vr1 - vn1; - vt2 = vr2 - vn2; - vt3 = vr3 - vn3; - - // relative rotational velocity - wr1 = (radi*omega[i][0] + radj*omega[j][0]); - wr2 = (radi*omega[i][1] + radj*omega[j][1]); - wr3 = (radi*omega[i][2] + radj*omega[j][2]); - - // relative tangential velocities - vtr1 = vt1 - (nz*wr2-ny*wr3); - vtr2 = vt2 - (nx*wr3-nz*wr1); - vtr3 = vt3 - (ny*wr1-nx*wr2); - vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; - vrel = sqrt(vrel); - - if (normal[itype][jtype] == JKR){ - F_pulloff = 3*M_PI*coh*Reff; - Fcrit = fabs(Fne + 2*F_pulloff); - } - else{ - Fcrit = fabs(Fne); - } - - //------------------------------ - //Tangential forces - //------------------------------ - k_tangential = tangential_coeffs[itype][jtype][0]; - if (normal[itype][jtype] != HOOKE){ - k_tangential *= a; - } - damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal; - - if (tangential_history){ - shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + - history[2]*history[2]); - - // tangential forces = history + tangential velocity damping - fs1 = -k_tangential*history[0] - damp_tangential*vtr1; - fs2 = -k_tangential*history[1] - damp_tangential*vtr2; - fs3 = -k_tangential*history[2] - damp_tangential*vtr3; - - // rescale frictional displacements and forces if needed - Fscrit = tangential_coeffs[itype][jtype][2] * Fcrit; - fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); - if (fs > Fscrit) { - if (shrmag != 0.0) { - history[0] = -1.0/k_tangential*(Fscrit*fs1/fs + damp_tangential*vtr1); - history[1] = -1.0/k_tangential*(Fscrit*fs2/fs + damp_tangential*vtr2); - history[2] = -1.0/k_tangential*(Fscrit*fs3/fs + damp_tangential*vtr3); - fs1 *= Fscrit/fs; - fs2 *= Fscrit/fs; - fs3 *= Fscrit/fs; - } else fs1 = fs2 = fs3 = 0.0; - } - } - else{ //Classic pair gran/hooke (no history) - fs = meff*damp_tangential*vrel; - if (vrel != 0.0) Ft = MIN(Fne,fs) / vrel; - else Ft = 0.0; - fs1 = -Ft*vtr1; - fs2 = -Ft*vtr2; - fs3 = -Ft*vtr3; - } - - //**************************************** - // Rolling resistance - //**************************************** - - if (roll[itype][jtype] != ROLL_NONE){ - relrot1 = omega[i][0] - omega[j][0]; - relrot2 = omega[i][1] - omega[j][1]; - relrot3 = omega[i][2] - omega[j][2]; - - // rolling velocity, see eq. 31 of Wang et al, Particuology v 23, p 49 (2015) - // This is different from the Marshall papers, which use the Bagi/Kuhn formulation - // for rolling velocity (see Wang et al for why the latter is wrong) - vrl1 = Reff*(relrot2*nz - relrot3*ny); //- 0.5*((radj-radi)/radsum)*vtr1; - vrl2 = Reff*(relrot3*nx - relrot1*nz); //- 0.5*((radj-radi)/radsum)*vtr2; - vrl3 = Reff*(relrot1*ny - relrot2*nx); //- 0.5*((radj-radi)/radsum)*vtr3; - vrlmag = sqrt(vrl1*vrl1+vrl2*vrl2+vrl3*vrl3); - if (vrlmag != 0.0) vrlmaginv = 1.0/vrlmag; - else vrlmaginv = 0.0; - - if (roll_history){ - int rhist0 = roll_history_index; - int rhist1 = rhist0 + 1; - int rhist2 = rhist1 + 1; - - // Rolling displacement - rollmag = sqrt(history[rhist0]*history[rhist0] + - history[rhist1]*history[rhist1] + - history[rhist2]*history[rhist2]); - - rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz; - - k_roll = roll_coeffs[itype][jtype][0]; - damp_roll = roll_coeffs[itype][jtype][1]; - fr1 = -k_roll*history[rhist0] - damp_roll*vrl1; - fr2 = -k_roll*history[rhist1] - damp_roll*vrl2; - fr3 = -k_roll*history[rhist2] - damp_roll*vrl3; - - // rescale frictional displacements and forces if needed - Frcrit = roll_coeffs[itype][jtype][2] * Fcrit; - - fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3); - if (fr > Frcrit) { - if (rollmag != 0.0) { - history[rhist0] = -1.0/k_roll*(Frcrit*fr1/fr + damp_roll*vrl1); - history[rhist1] = -1.0/k_roll*(Frcrit*fr2/fr + damp_roll*vrl2); - history[rhist2] = -1.0/k_roll*(Frcrit*fr3/fr + damp_roll*vrl3); - fr1 *= Frcrit/fr; - fr2 *= Frcrit/fr; - fr3 *= Frcrit/fr; - } else fr1 = fr2 = fr3 = 0.0; - } - } - else{ // - fr = meff*roll_coeffs[itype][jtype][1]*vrlmag; - if (vrlmag != 0.0) fr = MIN(Fne, fr) / vrlmag; - else fr = 0.0; - fr1 = -fr*vrl1; - fr2 = -fr*vrl2; - fr3 = -fr*vrl3; - } - } - - //**************************************** - // Twisting torque, including history effects - //**************************************** - if (twist[itype][jtype] != TWIST_NONE){ - magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of Marshall) - if (twist[itype][jtype] == TWIST_MARSHALL){ - k_twist = 0.5*k_tangential*a*a;; //eq 32 - damp_twist = 0.5*damp_tangential*a*a; - mu_twist = TWOTHIRDS*a; - } - else{ - k_twist = twist_coeffs[itype][jtype][0]; - damp_twist = twist_coeffs[itype][jtype][1]; - mu_twist = twist_coeffs[itype][jtype][2]; - } - if (twist_history){ - magtortwist = -k_twist*history[twist_history_index] - damp_twist*magtwist;//M_t torque (eq 30) - signtwist = (magtwist > 0) - (magtwist < 0); - Mtcrit = TWOTHIRDS*a*Fscrit;//critical torque (eq 44) - if (fabs(magtortwist) > Mtcrit) { - history[twist_history_index] = 1.0/k_twist*(Mtcrit*signtwist - damp_twist*magtwist); - magtortwist = -Mtcrit * signtwist; //eq 34 - } - } - else{ - if (magtwist > 0) magtortwist = -damp_twist*magtwist; - else magtortwist = 0; - } - } - - // set single_extra quantities - - svector[0] = fs1; - svector[1] = fs2; - svector[2] = fs3; - svector[3] = fs; - svector[4] = fr1; - svector[5] = fr2; - svector[6] = fr3; - svector[7] = fr; - svector[8] = magtortwist; - return 0.0; -} - -/* ---------------------------------------------------------------------- */ - -int PairGranularMulti::pack_forward_comm(int n, int *list, double *buf, - int pbc_flag, int *pbc) -{ - int i,j,m; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = mass_rigid[j]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -void PairGranularMulti::unpack_forward_comm(int n, int first, double *buf) -{ - int i,m,last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) - mass_rigid[i] = buf[m++]; -} - -/* ---------------------------------------------------------------------- - memory usage of local atom-based arrays - ------------------------------------------------------------------------- */ - -double PairGranularMulti::memory_usage() -{ - double bytes = nmax * sizeof(double); - return bytes; -} - -/* ---------------------------------------------------------------------- - mixing of Young's modulus (E) -------------------------------------------------------------------------- */ - -double PairGranularMulti::mix_stiffnessE(double Eii, double Ejj, double Gii, double Gjj) -{ - double poisii = Eii/(2.0*Gii) - 1.0; - double poisjj = Ejj/(2.0*Gjj) - 1.0; - return 1/((1-poisii*poisjj)/Eii+(1-poisjj*poisjj)/Ejj); -} - -/* ---------------------------------------------------------------------- - mixing of shear modulus (G) - ------------------------------------------------------------------------- */ - -double PairGranularMulti::mix_stiffnessG(double Eii, double Ejj, double Gii, double Gjj) -{ - double poisii = Eii/(2.0*Gii) - 1.0; - double poisjj = Ejj/(2.0*Gjj) - 1.0; - return 1/((2.0 -poisjj)/Gii+(2.0-poisjj)/Gjj); -} - -/* ---------------------------------------------------------------------- - mixing of everything else -------------------------------------------------------------------------- */ - -double PairGranularMulti::mix_geom(double valii, double valjj) -{ - return sqrt(valii*valjj); -} - - -/* ---------------------------------------------------------------------- - Compute pull-off distance (beyond contact) for a given radius and atom type -------------------------------------------------------------------------- */ - -double PairGranularMulti::pulloff_distance(double radius, int itype) -{ - double E, coh, a, delta_pulloff; - coh = normal_coeffs[itype][itype][3]; - E = mix_stiffnessE(normal_coeffs[itype][itype][0], normal_coeffs[itype][itype][0], - normal_coeffs[itype][itype][2], normal_coeffs[itype][itype][2]); - a = cbrt(9*M_PI*coh*radius*radius/(4*E)); - return a*a/radius - 2*sqrt(M_PI*coh*a/E); -} - diff --git a/src/GRANULAR/pair_granular_multi.h b/src/GRANULAR/pair_granular_multi.h deleted file mode 100644 index 95beb950f4..0000000000 --- a/src/GRANULAR/pair_granular_multi.h +++ /dev/null @@ -1,107 +0,0 @@ -/* ---------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifdef PAIR_CLASS - -PairStyle(granular/multi,PairGranularMulti) - -#else - -#ifndef LMP_PAIR_GRANULAR_MULTI_H -#define LMP_PAIR_GRANULAR_MULTI_H - -#include "pair.h" - -namespace LAMMPS_NS { - -class PairGranularMulti : public Pair { -public: - PairGranularMulti(class LAMMPS *); - virtual ~PairGranularMulti(); - virtual void compute(int, int); - virtual void settings(int, char **); - virtual void coeff(int, char **); - void init_style(); - double init_one(int, int); - void write_restart(FILE *); - void read_restart(FILE *); - void reset_dt(); - virtual double single(int, int, int, int, double, double, double, double &); - int pack_forward_comm(int, int *, double *, int, int *); - void unpack_forward_comm(int, int, double *); - double memory_usage(); - - protected: - double cut_global; - double dt; - int freeze_group_bit; - int use_history; - - int neighprev; - double *onerad_dynamic,*onerad_frozen; - double *maxrad_dynamic,*maxrad_frozen; - double **cut; - - class FixNeighHistory *fix_history; - - // storage of rigid body masses for use in granular interactions - - class Fix *fix_rigid; // ptr to rigid body fix, NULL if none - double *mass_rigid; // rigid mass for owned+ghost atoms - int nmax; // allocated size of mass_rigid - - virtual void allocate(); - -private: - int size_history; - - //Models - int **normal, **damping, **tangential, **roll, **twist; - - //History flags - int tangential_history, roll_history, twist_history; - - //Indices of history entries - int tangential_history_index; - int roll_history_index; - int twist_history_index; - - //Per-type coefficients, set in pair coeff command - double ***normal_coeffs; - double ***tangential_coeffs; - double ***roll_coeffs; - double ***twist_coeffs; - - //Optional user-specified global cutoff - double cutoff_global; - - double mix_stiffnessE(double Eii, double Ejj, double Gii, double Gjj); - double mix_stiffnessG(double Eii, double Ejj, double Gii, double Gjj); - double mix_geom(double valii, double valjj); - double pulloff_distance(double radius, int itype); -}; - -} - -#endif -#endif - -/* ERROR/WARNING messages: - -E: Illegal ... command - -Self-explanatory. Check the input script syntax and compare to the -documentation for the command. You can use -echo screen as a -command-line option when running LAMMPS to see the offending line. - - */