jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@15301 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2016-07-13 21:36:09 +00:00
parent b3aeb43a93
commit 6171c774aa
2 changed files with 106 additions and 46 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

141
src/GRANULAR/fix_wall_gran.cpp
View File

@ -44,7 +44,7 @@ enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY};
FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 10) error->all(FLERR,"Illegal fix wall/gran command");
if (narg < 11) error->all(FLERR,"Illegal fix wall/gran command");
if (!atom->sphere_flag)
error->all(FLERR,"Fix wall/gran requires atom style sphere");
@ -52,18 +52,28 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
restart_peratom = 1;
create_attribute = 1;
// wall/particle coefficients
// set interaction style
kn = force->numeric(FLERR,arg[3]);
if (strcmp(arg[4],"NULL") == 0) kt = kn * 2.0/7.0;
else kt = force->numeric(FLERR,arg[4]);
if (strcmp(arg[3],"hooke") == 0) pairstyle = HOOKE;
else if (strcmp(arg[3],"hooke/history") == 0) pairstyle = HOOKE_HISTORY;
else if (strcmp(arg[3],"hertz/history") == 0) pairstyle = HERTZ_HISTORY;
else error->all(FLERR,"Invalid fix wall/gran interaction style");
gamman = force->numeric(FLERR,arg[5]);
if (strcmp(arg[6],"NULL") == 0) gammat = 0.5 * gamman;
else gammat = force->numeric(FLERR,arg[6]);
history = 1;
if (pairstyle == HOOKE) history = 0;
xmu = force->numeric(FLERR,arg[7]);
int dampflag = force->inumeric(FLERR,arg[8]);
// particle/wall coefficients
kn = force->numeric(FLERR,arg[4]);
if (strcmp(arg[5],"NULL") == 0) kt = kn * 2.0/7.0;
else kt = force->numeric(FLERR,arg[5]);
gamman = force->numeric(FLERR,arg[6]);
if (strcmp(arg[7],"NULL") == 0) gammat = 0.5 * gamman;
else gammat = force->numeric(FLERR,arg[7]);
xmu = force->numeric(FLERR,arg[8]);
int dampflag = force->inumeric(FLERR,arg[9]);
if (dampflag == 0) gammat = 0.0;
if (kn < 0.0 || kt < 0.0 || gamman < 0.0 || gammat < 0.0 ||
@ -72,14 +82,14 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
// convert Kn and Kt from pressure units to force/distance^2 if Hertzian
if (force->pair_match("gran/hertz/history",1)) {
if (pairstyle == HERTZ_HISTORY) {
kn /= force->nktv2p;
kt /= force->nktv2p;
}
// wallstyle args
int iarg = 9;
int iarg = 10;
if (strcmp(arg[iarg],"xplane") == 0) {
if (narg < iarg+3) error->all(FLERR,"Illegal fix wall/gran command");
wallstyle = XPLANE;
@ -172,6 +182,9 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
atom->add_callback(0);
atom->add_callback(1);
nmax = 0;
mass_rigid = NULL;
// initialize as if particle is not touching wall
int nlocal = atom->nlocal;
@ -193,6 +206,7 @@ FixWallGran::~FixWallGran()
// delete locally stored arrays
memory->destroy(shear);
memory->destroy(mass_rigid);
}
/* ---------------------------------------------------------------------- */
@ -209,24 +223,19 @@ int FixWallGran::setmask()
void FixWallGran::init()
{
int i;
dt = update->dt;
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
// set pairstyle from granular pair style
// check for FixRigid so can extract rigid body masses
if (force->pair_match("gran/hooke",1))
pairstyle = HOOKE;
else if (force->pair_match("gran/hooke/history",1))
pairstyle = HOOKE_HISTORY;
else if (force->pair_match("gran/hooke/history/omp",1))
pairstyle = HOOKE_HISTORY;
else if (force->pair_match("gran/hertz/history",1))
pairstyle = HERTZ_HISTORY;
else if (force->pair_match("gran/hertz/history/omp",1))
pairstyle = HERTZ_HISTORY;
else error->all(FLERR,"Fix wall/gran is incompatible with Pair style");
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];
}
/* ---------------------------------------------------------------------- */
@ -246,11 +255,36 @@ void FixWallGran::setup(int vflag)
void FixWallGran::post_force(int vflag)
{
double vwall[3],dx,dy,dz,del1,del2,delxy,delr,rsq;
int i;
double dx,dy,dz,del1,del2,delxy,delr,rsq,meff;
double vwall[3];
// do not update shear history during setup
shearupdate = 1;
if (update->setupflag) shearupdate = 0;
// if just reneighbored:
// update rigid body masses for owned atoms if using FixRigid
// body[i] = which body atom I is in, -1 if none
// mass_body = mass of each rigid body
if (neighbor->ago == 0 && fix_rigid) {
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,"wall/gran: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;
}
}
// set position of wall to initial settings and velocity to 0.0
// if wiggle or shear, set wall position and velocity accordingly
@ -330,15 +364,24 @@ void FixWallGran::post_force(int vflag)
shear[i][2] = 0.0;
}
} else {
// meff = effective mass of sphere
// if I is part of rigid body, use body mass
meff = rmass[i];
if (fix_rigid && mass_rigid[i] > 0.0) meff = mass_rigid[i];
// inovke sphere/wall interaction
if (pairstyle == HOOKE)
hooke(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i],
radius[i],rmass[i]);
radius[i],meff);
else if (pairstyle == HOOKE_HISTORY)
hooke_history(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i],
radius[i],rmass[i],shear[i]);
radius[i],meff,shear[i]);
else if (pairstyle == HERTZ_HISTORY)
hertz_history(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i],
radius[i],rmass[i],shear[i]);
radius[i],meff,shear[i]);
}
}
}
@ -356,10 +399,10 @@ void FixWallGran::post_force_respa(int vflag, int ilevel, int iloop)
void FixWallGran::hooke(double rsq, double dx, double dy, double dz,
double *vwall, double *v,
double *f, double *omega, double *torque,
double radius, double mass)
double radius, double meff)
{
double r,vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
double wr1,wr2,wr3,meff,damp,ccel,vtr1,vtr2,vtr3,vrel;
double wr1,wr2,wr3,damp,ccel,vtr1,vtr2,vtr3,vrel;
double fn,fs,ft,fs1,fs2,fs3,fx,fy,fz,tor1,tor2,tor3,rinv,rsqinv;
r = sqrt(rsq);
@ -393,7 +436,6 @@ void FixWallGran::hooke(double rsq, double dx, double dy, double dz,
// normal forces = Hookian contact + normal velocity damping
meff = mass;
damp = meff*gamman*vnnr*rsqinv;
ccel = kn*(radius-r)*rinv - damp;
@ -441,10 +483,10 @@ void FixWallGran::hooke(double rsq, double dx, double dy, double dz,
void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz,
double *vwall, double *v,
double *f, double *omega, double *torque,
double radius, double mass, double *shear)
double radius, double meff, double *shear)
{
double r,vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
double wr1,wr2,wr3,meff,damp,ccel,vtr1,vtr2,vtr3,vrel;
double wr1,wr2,wr3,damp,ccel,vtr1,vtr2,vtr3,vrel;
double fn,fs,fs1,fs2,fs3,fx,fy,fz,tor1,tor2,tor3;
double shrmag,rsht,rinv,rsqinv;
@ -479,7 +521,6 @@ void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz,
// normal forces = Hookian contact + normal velocity damping
meff = mass;
damp = meff*gamman*vnnr*rsqinv;
ccel = kn*(radius-r)*rinv - damp;
@ -558,10 +599,10 @@ void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz,
void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz,
double *vwall, double *v,
double *f, double *omega, double *torque,
double radius, double mass, double *shear)
double radius, double meff, double *shear)
{
double r,vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
double wr1,wr2,wr3,meff,damp,ccel,vtr1,vtr2,vtr3,vrel;
double wr1,wr2,wr3,damp,ccel,vtr1,vtr2,vtr3,vrel;
double fn,fs,fs1,fs2,fs3,fx,fy,fz,tor1,tor2,tor3;
double shrmag,rsht,polyhertz,rinv,rsqinv;
@ -596,7 +637,6 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz,
// normal forces = Hertzian contact + normal velocity damping
meff = mass;
damp = meff*gamman*vnnr*rsqinv;
ccel = kn*(radius-r)*rinv - damp;
polyhertz = sqrt((radius-r)*radius);
@ -679,8 +719,9 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz,
double FixWallGran::memory_usage()
{
int nmax = atom->nmax;
double bytes = nmax * sizeof(int);
bytes += 3*nmax * sizeof(double);
double bytes = 0.0;
if (history) bytes += 3*nmax * sizeof(double); // shear history
if (fix_rigid) bytes += nmax * sizeof(int); // mass_rigid
return bytes;
}
@ -690,7 +731,7 @@ double FixWallGran::memory_usage()
void FixWallGran::grow_arrays(int nmax)
{
memory->grow(shear,nmax,3,"fix_wall_gran:shear");
if (history) memory->grow(shear,nmax,3,"fix_wall_gran:shear");
}
/* ----------------------------------------------------------------------
@ -699,9 +740,11 @@ void FixWallGran::grow_arrays(int nmax)
void FixWallGran::copy_arrays(int i, int j, int delflag)
{
shear[j][0] = shear[i][0];
shear[j][1] = shear[i][1];
shear[j][2] = shear[i][2];
if (history) {
shear[j][0] = shear[i][0];
shear[j][1] = shear[i][1];
shear[j][2] = shear[i][2];
}
}
/* ----------------------------------------------------------------------
@ -710,7 +753,7 @@ void FixWallGran::copy_arrays(int i, int j, int delflag)
void FixWallGran::set_arrays(int i)
{
shear[i][0] = shear[i][1] = shear[i][2] = 0.0;
if (history) shear[i][0] = shear[i][1] = shear[i][2] = 0.0;
}
/* ----------------------------------------------------------------------
@ -719,6 +762,8 @@ void FixWallGran::set_arrays(int i)
int FixWallGran::pack_exchange(int i, double *buf)
{
if (!history) return 0;
buf[0] = shear[i][0];
buf[1] = shear[i][1];
buf[2] = shear[i][2];
@ -731,6 +776,8 @@ int FixWallGran::pack_exchange(int i, double *buf)
int FixWallGran::unpack_exchange(int nlocal, double *buf)
{
if (!history) return 0;
shear[nlocal][0] = buf[0];
shear[nlocal][1] = buf[1];
shear[nlocal][2] = buf[2];
@ -743,6 +790,8 @@ int FixWallGran::unpack_exchange(int nlocal, double *buf)
int FixWallGran::pack_restart(int i, double *buf)
{
if (!history) return 0;
int m = 0;
buf[m++] = 4;
buf[m++] = shear[i][0];
@ -759,6 +808,8 @@ void FixWallGran::unpack_restart(int nlocal, int nth)
{
double **extra = atom->extra;
if (!history) return;
// skip to Nth set of extra values
int m = 0;
@ -776,6 +827,7 @@ void FixWallGran::unpack_restart(int nlocal, int nth)
int FixWallGran::maxsize_restart()
{
if (!history) return 0;
return 4;
}
@ -785,6 +837,7 @@ int FixWallGran::maxsize_restart()
int FixWallGran::size_restart(int nlocal)
{
if (!history) return 0;
return 4;
}

11
src/GRANULAR/fix_wall_gran.h
View File

@ -47,7 +47,7 @@ class FixWallGran : public Fix {
void reset_dt();
protected:
int wallstyle,pairstyle,wiggle,wshear,axis;
int wallstyle,pairstyle,history,wiggle,wshear,axis;
double kn,kt,gamman,gammat,xmu;
double lo,hi,cylradius;
double amplitude,period,omega,vshear;
@ -55,10 +55,17 @@ class FixWallGran : public Fix {
int nlevels_respa;
int time_origin;
int *touch;
// shear history values
double **shear;
int shearupdate;
// 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
void hooke(double, double, double, double, double *,
double *, double *, double *, double *, double, double);
void hooke_history(double, double, double, double, double *,