From 7e0d44e0ca471495e5a8dd4bb90040fd14a6aa84 Mon Sep 17 00:00:00 2001 From: Joel Clemmer Date: Tue, 30 Mar 2021 20:51:23 -0600 Subject: [PATCH 01/15] Adding no_attraction flag to granular pair and wall styles --- doc/src/fix_wall_gran.rst | 10 ++++++++-- doc/src/fix_wall_gran_region.rst | 10 +++++++++- doc/src/pair_gran.rst | 15 ++++++++++++++- doc/src/pair_granular.rst | 6 ++++++ src/GRANULAR/fix_wall_gran.cpp | 14 +++++++++++++- src/GRANULAR/fix_wall_gran.h | 1 + src/GRANULAR/pair_gran_hertz_history.cpp | 2 ++ src/GRANULAR/pair_gran_hooke.cpp | 2 ++ src/GRANULAR/pair_gran_hooke_history.cpp | 12 +++++++++++- src/GRANULAR/pair_gran_hooke_history.h | 1 + src/GRANULAR/pair_granular.cpp | 14 +++++++++++++- src/GRANULAR/pair_granular.h | 1 + src/KOKKOS/pair_gran_hooke_history_kokkos.cpp | 1 + 13 files changed, 82 insertions(+), 7 deletions(-) diff --git a/doc/src/fix_wall_gran.rst b/doc/src/fix_wall_gran.rst index 1eefc4aeae..73fa7e367b 100644 --- a/doc/src/fix_wall_gran.rst +++ b/doc/src/fix_wall_gran.rst @@ -45,7 +45,7 @@ Syntax radius = cylinder radius (distance units) * zero or more keyword/value pairs may be appended to args -* keyword = *wiggle* or *shear* or *contacts* +* keyword = *wiggle* or *shear* or *contacts* or *no_attraction* .. parsed-literal:: @@ -58,6 +58,8 @@ Syntax vshear = magnitude of shear velocity (velocity units) *contacts* value = none generate contact information for each particle + *no_attraction* value = none + turn off possibility for attractive interactions Examples @@ -175,8 +177,12 @@ the clockwise direction for *vshear* > 0 or counter-clockwise for *vshear* < 0. In this case, *vshear* is the tangential velocity of the wall at whatever *radius* has been defined. +If a particle is moving away from the wall while in contact, there +is a possibility that the particle could experience an effective attractive +force due to damping. If the *no_attraction* keyword is used, this fix +will zero out the normal component of the force if there is an effective +attractive force. This keyword cannot be used with the JKR or DMT models. -Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" This fix writes the shear friction state of atoms interacting with the diff --git a/doc/src/fix_wall_gran_region.rst b/doc/src/fix_wall_gran_region.rst index 9ed5bb7bcc..356c8427b0 100644 --- a/doc/src/fix_wall_gran_region.rst +++ b/doc/src/fix_wall_gran_region.rst @@ -36,12 +36,14 @@ Syntax * wallstyle = region (see :doc:`fix wall/gran ` for options for other kinds of walls) * region-ID = region whose boundary will act as wall -* keyword = *contacts* +* keyword = *contacts* or *no_attraction* .. parsed-literal:: *contacts* value = none generate contact information for each particle + *no_attraction* value = none + turn off possibility for attractive interactions Examples """""""" @@ -199,6 +201,12 @@ values for the 6 wall/particle coefficients than for particle/particle interactions. E.g. if you wish to model the wall as a different material. +If a particle is moving away from the wall while in contact, there +is a possibility that the particle could experience an effective attractive +force due to damping. If the *no_attraction* keyword is used, this fix +will zero out the normal component of the force if there is an effective +attractive force. This keyword cannot be used with the JKR or DMT models. + Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" diff --git a/doc/src/pair_gran.rst b/doc/src/pair_gran.rst index d3f87839e8..aa2e42d339 100644 --- a/doc/src/pair_gran.rst +++ b/doc/src/pair_gran.rst @@ -26,7 +26,7 @@ Syntax .. code-block:: LAMMPS - pair_style style Kn Kt gamma_n gamma_t xmu dampflag + pair_style style Kn Kt gamma_n gamma_t xmu dampflag keyword * style = *gran/hooke* or *gran/hooke/history* or *gran/hertz/history* * Kn = elastic constant for normal particle repulsion (force/distance units or pressure units - see discussion below) @@ -36,6 +36,13 @@ Syntax * xmu = static yield criterion (unitless value between 0.0 and 1.0e4) * dampflag = 0 or 1 if tangential damping force is excluded or included +* keyword = *no_attraction* + + .. parsed-literal:: + + *no_attraction* value = none + turn off possibility for attractive interactions + .. note:: Versions of LAMMPS before 9Jan09 had different style names for @@ -208,6 +215,12 @@ potential is used as a sub-style of :doc:`pair_style hybrid `, then pair_coeff command to determine which atoms interact via a granular potential. +If two particles are moving away from each other while in contact, there +is a possibility that the particles could experience an effective attractive +force due to damping. If the *no_attraction* keyword is used, this fix +will zero out the normal component of the force if there is an effective +attractive force. + ---------- .. include:: accel_styles.rst diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index 80663b7d49..cb65a58777 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -657,6 +657,12 @@ then LAMMPS will use that cutoff for the specified atom type combination, and automatically set pairwise cutoffs for the remaining atom types. +If two particles are moving away from each other while in contact, there +is a possibility that the particles could experience an effective attractive +force due to damping. If the *no_attraction* keyword is used, this fix +will zero out the normal component of the force if there is an effective +attractive force. This keyword cannot be used with the JKR or DMT models. + ---------- .. include:: accel_styles.rst diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index d18b72aa65..9c17a323c5 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -347,6 +347,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : wiggle = 0; wshear = 0; peratom_flag = 0; + noattraction_flag = 0; while (iarg < narg) { if (strcmp(arg[iarg],"wiggle") == 0) { @@ -373,6 +374,13 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : size_peratom_cols = 8; peratom_freq = 1; iarg += 1; + } else if (strcmp(arg[iarg],"no_attraction") == 0) { + noattraction_flag = 1; + if(normal_model == JKR) + error->all(FLERR,"Cannot turn off attraction with JRK model"); + if(normal_model == DMT) + error->all(FLERR,"Cannot turn off attraction with DMT model"); + iarg += 1; } else error->all(FLERR,"Illegal fix wall/gran command"); } @@ -761,6 +769,7 @@ void FixWallGran::hooke(double rsq, double dx, double dy, double dz, damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radius-r)*rinv - damp; + if(noattraction_flag and ccel < 0.0) ccel = 0.0; // relative velocities @@ -853,6 +862,7 @@ void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz, damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radius-r)*rinv - damp; + if(noattraction_flag and ccel < 0.0) ccel = 0.0; // relative velocities @@ -984,7 +994,8 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz, if (rwall == 0.0) polyhertz = sqrt((radius-r)*radius); else polyhertz = sqrt((radius-r)*radius*rwall/(rwall+radius)); ccel *= polyhertz; - + if(noattraction_flag and ccel < 0.0) ccel = 0.0; + // relative velocities vtr1 = vt1 - (dz*wr2-dy*wr3); @@ -1177,6 +1188,7 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz, Fdamp = -damp_normal_prefactor*vnnr; Fntot = Fne + Fdamp; + if(noattraction_flag and Fntot < 0.0) Fntot = 0.0; //**************************************** // tangential force, including history effects diff --git a/src/GRANULAR/fix_wall_gran.h b/src/GRANULAR/fix_wall_gran.h index a81a4c787c..3963b22a2c 100644 --- a/src/GRANULAR/fix_wall_gran.h +++ b/src/GRANULAR/fix_wall_gran.h @@ -69,6 +69,7 @@ class FixWallGran : public Fix { // for granular model choices int normal_model, damping_model; int tangential_model, roll_model, twist_model; + int noattraction_flag; // history flags int normal_history, tangential_history, roll_history, twist_history; diff --git a/src/GRANULAR/pair_gran_hertz_history.cpp b/src/GRANULAR/pair_gran_hertz_history.cpp index e7193140b9..7da0c6142d 100644 --- a/src/GRANULAR/pair_gran_hertz_history.cpp +++ b/src/GRANULAR/pair_gran_hertz_history.cpp @@ -181,6 +181,7 @@ void PairGranHertzHistory::compute(int eflag, int vflag) ccel = kn*(radsum-r)*rinv - damp; polyhertz = sqrt((radsum-r)*radi*radj / radsum); ccel *= polyhertz; + if(noattractionflag and ccel < 0.0) ccel = 0.0; // relative velocities @@ -388,6 +389,7 @@ double PairGranHertzHistory::single(int i, int j, int /*itype*/, int /*jtype*/, ccel = kn*(radsum-r)*rinv - damp; polyhertz = sqrt((radsum-r)*radi*radj / radsum); ccel *= polyhertz; + if(noattractionflag and ccel < 0.0) ccel = 0.0; // relative velocities diff --git a/src/GRANULAR/pair_gran_hooke.cpp b/src/GRANULAR/pair_gran_hooke.cpp index 3875d12a0f..76946701cb 100644 --- a/src/GRANULAR/pair_gran_hooke.cpp +++ b/src/GRANULAR/pair_gran_hooke.cpp @@ -158,6 +158,7 @@ void PairGranHooke::compute(int eflag, int vflag) damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; + if(noattractionflag and ccel < 0.0) ccel = 0.0; // relative velocities @@ -299,6 +300,7 @@ double PairGranHooke::single(int i, int j, int /*itype*/, int /*jtype*/, double damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; + if(noattractionflag and ccel < 0.0) ccel = 0.0; // relative velocities diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp index 4deb828d76..4893be4084 100644 --- a/src/GRANULAR/pair_gran_hooke_history.cpp +++ b/src/GRANULAR/pair_gran_hooke_history.cpp @@ -237,6 +237,7 @@ void PairGranHookeHistory::compute(int eflag, int vflag) damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; + if(noattractionflag and ccel < 0.0) ccel = 0.0; // relative velocities @@ -356,7 +357,7 @@ void PairGranHookeHistory::allocate() void PairGranHookeHistory::settings(int narg, char **arg) { - if (narg != 6) error->all(FLERR,"Illegal pair_style command"); + if (narg != 6 and narg != 7) error->all(FLERR,"Illegal pair_style command"); kn = utils::numeric(FLERR,arg[0],false,lmp); if (strcmp(arg[1],"NULL") == 0) kt = kn * 2.0/7.0; @@ -369,6 +370,14 @@ void PairGranHookeHistory::settings(int narg, char **arg) xmu = utils::numeric(FLERR,arg[4],false,lmp); dampflag = utils::inumeric(FLERR,arg[5],false,lmp); if (dampflag == 0) gammat = 0.0; + + noattractionflag = 0; + if(narg == 7) { + if(strcmp(arg[6], "no_attraction")) + noattractionflag = 1; + else + error->all(FLERR,"Illegal pair_style command"); + } if (kn < 0.0 || kt < 0.0 || gamman < 0.0 || gammat < 0.0 || xmu < 0.0 || xmu > 10000.0 || dampflag < 0 || dampflag > 1) @@ -686,6 +695,7 @@ double PairGranHookeHistory::single(int i, int j, int /*itype*/, int /*jtype*/, damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; + if(noattractionflag and ccel < 0.0) ccel = 0.0; // relative velocities diff --git a/src/GRANULAR/pair_gran_hooke_history.h b/src/GRANULAR/pair_gran_hooke_history.h index c27ce8a9af..bf40097f9d 100644 --- a/src/GRANULAR/pair_gran_hooke_history.h +++ b/src/GRANULAR/pair_gran_hooke_history.h @@ -49,6 +49,7 @@ class PairGranHookeHistory : public Pair { double dt; int freeze_group_bit; int history; + int noattractionflag; int neighprev; double *onerad_dynamic,*onerad_frozen; diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index cc210138eb..7c1a1f2b29 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -793,6 +793,7 @@ void PairGranular::coeff(int narg, char **arg) roll_model_one = ROLL_NONE; twist_model_one = TWIST_NONE; damping_model_one = VISCOELASTIC; + noattractionflag = 0; int iarg = 2; while (iarg < narg) { @@ -1030,7 +1031,18 @@ void PairGranular::coeff(int narg, char **arg) count++; } } - + + if(noattraction flag) { + for(int i = 1; i <= ntypes; i++) { + for(int j = 1; j <= ntypes; j++) { + if(normal_model[i][j] == JKR) + error->all(FLERR,"Cannot turn off attraction with JKR pairstyle"); + if(normal_model[i][j] == DMT) + error->all(FLERR,"Cannot turn off attraction with DMT pairstyle"); + } + } + } + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } diff --git a/src/GRANULAR/pair_granular.h b/src/GRANULAR/pair_granular.h index 7ef4f2af98..1e4e8c8fc4 100644 --- a/src/GRANULAR/pair_granular.h +++ b/src/GRANULAR/pair_granular.h @@ -70,6 +70,7 @@ class PairGranular : public Pair { // model choices int **normal_model, **damping_model; int **tangential_model, **roll_model, **twist_model; + int **noattractionflag; // history flags int normal_history, tangential_history, roll_history, twist_history; diff --git a/src/KOKKOS/pair_gran_hooke_history_kokkos.cpp b/src/KOKKOS/pair_gran_hooke_history_kokkos.cpp index 81c82ab826..caa19ffdef 100644 --- a/src/KOKKOS/pair_gran_hooke_history_kokkos.cpp +++ b/src/KOKKOS/pair_gran_hooke_history_kokkos.cpp @@ -392,6 +392,7 @@ void PairGranHookeHistoryKokkos::operator()(TagPairGranHookeHistoryC F_FLOAT damp = meff*gamman*vnnr*rsqinv; F_FLOAT ccel = kn*(radsum-r)*rinv - damp; + if(noattractionflag & ccel < 0.0) ccel = 0.0; // relative velocities From 2988aa2d1f6cc3d6a67bc7b03021d898656a30f9 Mon Sep 17 00:00:00 2001 From: Dan Bolintineanu Date: Tue, 30 Mar 2021 21:55:06 -0600 Subject: [PATCH 02/15] Fixes bug where accumulated tangential and rolling displacements are not rotated correctly. Pointed out by Deng Pan on LAMMPS mailing list 3/26/21 --- src/GRANULAR/pair_granular.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index cc210138eb..bcf262aa2c 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -446,9 +446,9 @@ void PairGranular::compute(int eflag, int vflag) if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_FORCE || tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE_FORCE) - frameupdate = fabs(rsht) < EPSILON*Fscrit; + frameupdate = fabs(rsht) > EPSILON*Fscrit; else - frameupdate = fabs(rsht)*k_tangential < EPSILON*Fscrit; + frameupdate = fabs(rsht)*k_tangential > EPSILON*Fscrit; if (frameupdate) { shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + history[2]*history[2]); @@ -568,7 +568,7 @@ void PairGranular::compute(int eflag, int vflag) if (historyupdate) { rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz; - frameupdate = fabs(rolldotn)*k_roll < EPSILON*Frcrit; + frameupdate = fabs(rolldotn)*k_roll > EPSILON*Frcrit; if (frameupdate) { // rotate into tangential plane rollmag = sqrt(history[rhist0]*history[rhist0] + history[rhist1]*history[rhist1] + From d36894ccf5691ca9ba6e85d264a3cb435526ce44 Mon Sep 17 00:00:00 2001 From: Joel Clemmer Date: Tue, 30 Mar 2021 22:16:47 -0600 Subject: [PATCH 03/15] Adding pertype flag to pair granular --- doc/src/fix_wall_gran.rst | 2 +- doc/src/fix_wall_gran_region.rst | 2 +- doc/src/pair_gran.rst | 2 +- doc/src/pair_granular.rst | 8 ++++++++ src/GRANULAR/pair_granular.cpp | 34 ++++++++++++++++++++------------ 5 files changed, 32 insertions(+), 16 deletions(-) diff --git a/doc/src/fix_wall_gran.rst b/doc/src/fix_wall_gran.rst index 73fa7e367b..bfe06bf7a0 100644 --- a/doc/src/fix_wall_gran.rst +++ b/doc/src/fix_wall_gran.rst @@ -59,7 +59,7 @@ Syntax *contacts* value = none generate contact information for each particle *no_attraction* value = none - turn off possibility for attractive interactions + turn off possibility of attractive interactions Examples diff --git a/doc/src/fix_wall_gran_region.rst b/doc/src/fix_wall_gran_region.rst index 356c8427b0..00eb19ba1c 100644 --- a/doc/src/fix_wall_gran_region.rst +++ b/doc/src/fix_wall_gran_region.rst @@ -43,7 +43,7 @@ Syntax *contacts* value = none generate contact information for each particle *no_attraction* value = none - turn off possibility for attractive interactions + turn off possibility of attractive interactions Examples """""""" diff --git a/doc/src/pair_gran.rst b/doc/src/pair_gran.rst index aa2e42d339..8145f7fb7a 100644 --- a/doc/src/pair_gran.rst +++ b/doc/src/pair_gran.rst @@ -41,7 +41,7 @@ Syntax .. parsed-literal:: *no_attraction* value = none - turn off possibility for attractive interactions + turn off possibility of attractive interactions .. note:: diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index cb65a58777..0076994e25 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -623,6 +623,14 @@ Finally, the twisting torque on each particle is given by: ---------- +If two particles are moving away from each other while in contact, there +is a possibility that the particles could experience an effective attractive +force due to damping. If the optional *no_attraction* keyword is used, this fix +will zero out the normal component of the force if there is an effective +attractive force. This keyword cannot be used with the JKR or DMT models. + +---------- + The *granular* pair style can reproduce the behavior of the *pair gran/\** styles with the appropriate settings (some very minor differences can be expected due to corrections in diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 7c1a1f2b29..256acd5a6a 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -130,6 +130,7 @@ PairGranular::~PairGranular() memory->destroy(tangential_model); memory->destroy(roll_model); memory->destroy(twist_model); + memory->destroy(noattraction_flag); delete [] onerad_dynamic; delete [] onerad_frozen; @@ -370,6 +371,7 @@ void PairGranular::compute(int eflag, int vflag) Fdamp = -damp_normal_prefactor*vnnr; Fntot = Fne + Fdamp; + if(noattraction_flag[itype][jtype] and Fntot < 0.0) Fntot = 0.0; //**************************************** // tangential force, including history effects @@ -740,6 +742,7 @@ void PairGranular::allocate() memory->create(tangential_model,n+1,n+1,"pair:tangential_model"); memory->create(roll_model,n+1,n+1,"pair:roll_model"); memory->create(twist_model,n+1,n+1,"pair:twist_model"); + memory->create(noattraction_flag,n+1,n+1,"pair:noattraction_flag"); onerad_dynamic = new double[n+1]; onerad_frozen = new double[n+1]; @@ -793,8 +796,8 @@ void PairGranular::coeff(int narg, char **arg) roll_model_one = ROLL_NONE; twist_model_one = TWIST_NONE; damping_model_one = VISCOELASTIC; - noattractionflag = 0; - + int na_flag = 0; + int iarg = 2; while (iarg < narg) { if (strcmp(arg[iarg], "hooke") == 0) { @@ -968,6 +971,9 @@ void PairGranular::coeff(int narg, char **arg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); cutoff_one = utils::numeric(FLERR,arg[iarg+1],false,lmp); iarg += 2; + } else if (strcmp(arg[iarg], "no_attraction") == 0) { + na_flag = 1; + iarg += 1; } else error->all(FLERR, "Illegal pair coeff command"); } @@ -985,6 +991,12 @@ void PairGranular::coeff(int narg, char **arg) 27.467*powint(cor,4)-18.022*powint(cor,5)+4.8218*powint(cor,6); } else damp = normal_coeffs_one[1]; + if(na_flag and normal_model_one == JKR) + error->all(FLERR,"Cannot turn off attraction with JKR pairstyle"); + + if(na_flag and normal_model_one == DMT) + error->all(FLERR,"Cannot turn off attraction with DMT pairstyle"); + for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { normal_model[i][j] = normal_model[j][i] = normal_model_one; @@ -1027,21 +1039,13 @@ void PairGranular::coeff(int narg, char **arg) cutoff_type[i][j] = cutoff_type[j][i] = cutoff_one; + noattraction_flag[i][j] = na_flag; + setflag[i][j] = 1; count++; } } - - if(noattraction flag) { - for(int i = 1; i <= ntypes; i++) { - for(int j = 1; j <= ntypes; j++) { - if(normal_model[i][j] == JKR) - error->all(FLERR,"Cannot turn off attraction with JKR pairstyle"); - if(normal_model[i][j] == DMT) - error->all(FLERR,"Cannot turn off attraction with DMT pairstyle"); - } - } - } + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } @@ -1315,6 +1319,7 @@ void PairGranular::write_restart(FILE *fp) fwrite(&tangential_model[i][j],sizeof(int),1,fp); fwrite(&roll_model[i][j],sizeof(int),1,fp); fwrite(&twist_model[i][j],sizeof(int),1,fp); + fwrite(&noattraction_flag[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); @@ -1345,6 +1350,7 @@ void PairGranular::read_restart(FILE *fp) utils::sfread(FLERR,&tangential_model[i][j],sizeof(int),1,fp,nullptr,error); utils::sfread(FLERR,&roll_model[i][j],sizeof(int),1,fp,nullptr,error); utils::sfread(FLERR,&twist_model[i][j],sizeof(int),1,fp,nullptr,error); + utils::sfread(FLERR,&noattraction_flag[i][j],sizeof(int),1,fp,nullptr,error); utils::sfread(FLERR,normal_coeffs[i][j],sizeof(double),4,fp,nullptr,error); utils::sfread(FLERR,tangential_coeffs[i][j],sizeof(double),3,fp,nullptr,error); utils::sfread(FLERR,roll_coeffs[i][j],sizeof(double),3,fp,nullptr,error); @@ -1356,6 +1362,7 @@ void PairGranular::read_restart(FILE *fp) MPI_Bcast(&tangential_model[i][j],1,MPI_INT,0,world); MPI_Bcast(&roll_model[i][j],1,MPI_INT,0,world); MPI_Bcast(&twist_model[i][j],1,MPI_INT,0,world); + MPI_Bcast(&noattraction_flag[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); @@ -1541,6 +1548,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, Fdamp = -damp_normal_prefactor*vnnr; Fntot = Fne + Fdamp; + if(noattraction_flag[itype][jtype] and Fntot < 0.0) Fntot = 0.0; jnum = list->numneigh[i]; jlist = list->firstneigh[i]; From 8b2676a103c8fc220fe404c547d1e40ea3d78d5f Mon Sep 17 00:00:00 2001 From: Joel Clemmer Date: Wed, 31 Mar 2021 16:57:13 -0600 Subject: [PATCH 04/15] Moving arguments in gran wall, renaming toarg limit_damping --- src/GRANULAR/fix_wall_gran.cpp | 31 ++++++++++++------- src/GRANULAR/fix_wall_gran.h | 2 +- src/GRANULAR/pair_gran_hertz_history.cpp | 4 +-- src/GRANULAR/pair_gran_hooke.cpp | 4 +-- src/GRANULAR/pair_gran_hooke_history.cpp | 10 +++--- src/GRANULAR/pair_gran_hooke_history.h | 2 +- src/GRANULAR/pair_granular.cpp | 30 +++++++++--------- src/GRANULAR/pair_granular.h | 2 +- src/KOKKOS/pair_gran_hooke_history_kokkos.cpp | 2 +- 9 files changed, 47 insertions(+), 40 deletions(-) diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index 9c17a323c5..c9d00ccdd3 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -116,6 +116,12 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : kt /= force->nktv2p; } iarg = 10; + + if (strcmp(arg[iarg],"limit_damping") == 0) { + limit_damping = 1; + iarg += 1; + } + } else { iarg = 4; damping_model = VISCOELASTIC; @@ -292,6 +298,9 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : strcmp(arg[iarg], "zcylinder") == 0 || strcmp(arg[iarg], "region") == 0) { break; + } else if (strcmp(arg[iarg],"limit_damping") == 0) { + limit_damping = 1; + iarg += 1; } else { error->all(FLERR, "Illegal fix wall/gran command"); } @@ -301,6 +310,11 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : if (tangential_model == TANGENTIAL_MINDLIN_RESCALE) size_history += 1; } + if(normal_model == JKR) + error->all(FLERR,"Cannot limit damping with JRK model"); + if(normal_model == DMT) + error->all(FLERR,"Cannot limit damping with DMT model"); + // wallstyle args idregion = nullptr; @@ -347,7 +361,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : wiggle = 0; wshear = 0; peratom_flag = 0; - noattraction_flag = 0; + limit_damping = 0; while (iarg < narg) { if (strcmp(arg[iarg],"wiggle") == 0) { @@ -374,13 +388,6 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : size_peratom_cols = 8; peratom_freq = 1; iarg += 1; - } else if (strcmp(arg[iarg],"no_attraction") == 0) { - noattraction_flag = 1; - if(normal_model == JKR) - error->all(FLERR,"Cannot turn off attraction with JRK model"); - if(normal_model == DMT) - error->all(FLERR,"Cannot turn off attraction with DMT model"); - iarg += 1; } else error->all(FLERR,"Illegal fix wall/gran command"); } @@ -769,7 +776,7 @@ void FixWallGran::hooke(double rsq, double dx, double dy, double dz, damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radius-r)*rinv - damp; - if(noattraction_flag and ccel < 0.0) ccel = 0.0; + if(limit_damping and ccel < 0.0) ccel = 0.0; // relative velocities @@ -862,7 +869,7 @@ void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz, damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radius-r)*rinv - damp; - if(noattraction_flag and ccel < 0.0) ccel = 0.0; + if(limit_damping and ccel < 0.0) ccel = 0.0; // relative velocities @@ -994,7 +1001,7 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz, if (rwall == 0.0) polyhertz = sqrt((radius-r)*radius); else polyhertz = sqrt((radius-r)*radius*rwall/(rwall+radius)); ccel *= polyhertz; - if(noattraction_flag and ccel < 0.0) ccel = 0.0; + if(limit_damping and ccel < 0.0) ccel = 0.0; // relative velocities @@ -1188,7 +1195,7 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz, Fdamp = -damp_normal_prefactor*vnnr; Fntot = Fne + Fdamp; - if(noattraction_flag and Fntot < 0.0) Fntot = 0.0; + if(limit_damping and Fntot < 0.0) Fntot = 0.0; //**************************************** // tangential force, including history effects diff --git a/src/GRANULAR/fix_wall_gran.h b/src/GRANULAR/fix_wall_gran.h index 3963b22a2c..473f26cca4 100644 --- a/src/GRANULAR/fix_wall_gran.h +++ b/src/GRANULAR/fix_wall_gran.h @@ -69,7 +69,7 @@ class FixWallGran : public Fix { // for granular model choices int normal_model, damping_model; int tangential_model, roll_model, twist_model; - int noattraction_flag; + int limit_damping; // history flags int normal_history, tangential_history, roll_history, twist_history; diff --git a/src/GRANULAR/pair_gran_hertz_history.cpp b/src/GRANULAR/pair_gran_hertz_history.cpp index 7da0c6142d..40c6c3a41d 100644 --- a/src/GRANULAR/pair_gran_hertz_history.cpp +++ b/src/GRANULAR/pair_gran_hertz_history.cpp @@ -181,7 +181,7 @@ void PairGranHertzHistory::compute(int eflag, int vflag) ccel = kn*(radsum-r)*rinv - damp; polyhertz = sqrt((radsum-r)*radi*radj / radsum); ccel *= polyhertz; - if(noattractionflag and ccel < 0.0) ccel = 0.0; + if(limit_damping and ccel < 0.0) ccel = 0.0; // relative velocities @@ -389,7 +389,7 @@ double PairGranHertzHistory::single(int i, int j, int /*itype*/, int /*jtype*/, ccel = kn*(radsum-r)*rinv - damp; polyhertz = sqrt((radsum-r)*radi*radj / radsum); ccel *= polyhertz; - if(noattractionflag and ccel < 0.0) ccel = 0.0; + if(limit_damping and ccel < 0.0) ccel = 0.0; // relative velocities diff --git a/src/GRANULAR/pair_gran_hooke.cpp b/src/GRANULAR/pair_gran_hooke.cpp index 76946701cb..0e3ffc7293 100644 --- a/src/GRANULAR/pair_gran_hooke.cpp +++ b/src/GRANULAR/pair_gran_hooke.cpp @@ -158,7 +158,7 @@ void PairGranHooke::compute(int eflag, int vflag) damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; - if(noattractionflag and ccel < 0.0) ccel = 0.0; + if(limit_damping and ccel < 0.0) ccel = 0.0; // relative velocities @@ -300,7 +300,7 @@ double PairGranHooke::single(int i, int j, int /*itype*/, int /*jtype*/, double damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; - if(noattractionflag and ccel < 0.0) ccel = 0.0; + if(limit_damping and ccel < 0.0) ccel = 0.0; // relative velocities diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp index 4893be4084..d3c1cf7672 100644 --- a/src/GRANULAR/pair_gran_hooke_history.cpp +++ b/src/GRANULAR/pair_gran_hooke_history.cpp @@ -237,7 +237,7 @@ void PairGranHookeHistory::compute(int eflag, int vflag) damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; - if(noattractionflag and ccel < 0.0) ccel = 0.0; + if(limit_damping and ccel < 0.0) ccel = 0.0; // relative velocities @@ -371,10 +371,10 @@ void PairGranHookeHistory::settings(int narg, char **arg) dampflag = utils::inumeric(FLERR,arg[5],false,lmp); if (dampflag == 0) gammat = 0.0; - noattractionflag = 0; + limit_damping = 0; if(narg == 7) { - if(strcmp(arg[6], "no_attraction")) - noattractionflag = 1; + if(strcmp(arg[6], "limit_damping")) + limit_damping = 1; else error->all(FLERR,"Illegal pair_style command"); } @@ -695,7 +695,7 @@ double PairGranHookeHistory::single(int i, int j, int /*itype*/, int /*jtype*/, damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; - if(noattractionflag and ccel < 0.0) ccel = 0.0; + if(limit_damping and ccel < 0.0) ccel = 0.0; // relative velocities diff --git a/src/GRANULAR/pair_gran_hooke_history.h b/src/GRANULAR/pair_gran_hooke_history.h index bf40097f9d..f511d30237 100644 --- a/src/GRANULAR/pair_gran_hooke_history.h +++ b/src/GRANULAR/pair_gran_hooke_history.h @@ -49,7 +49,7 @@ class PairGranHookeHistory : public Pair { double dt; int freeze_group_bit; int history; - int noattractionflag; + int limit_damping; int neighprev; double *onerad_dynamic,*onerad_frozen; diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 256acd5a6a..dd22f03fa8 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -130,7 +130,7 @@ PairGranular::~PairGranular() memory->destroy(tangential_model); memory->destroy(roll_model); memory->destroy(twist_model); - memory->destroy(noattraction_flag); + memory->destroy(limit_damping); delete [] onerad_dynamic; delete [] onerad_frozen; @@ -371,7 +371,7 @@ void PairGranular::compute(int eflag, int vflag) Fdamp = -damp_normal_prefactor*vnnr; Fntot = Fne + Fdamp; - if(noattraction_flag[itype][jtype] and Fntot < 0.0) Fntot = 0.0; + if(limit_damping[itype][jtype] and Fntot < 0.0) Fntot = 0.0; //**************************************** // tangential force, including history effects @@ -742,7 +742,7 @@ void PairGranular::allocate() memory->create(tangential_model,n+1,n+1,"pair:tangential_model"); memory->create(roll_model,n+1,n+1,"pair:roll_model"); memory->create(twist_model,n+1,n+1,"pair:twist_model"); - memory->create(noattraction_flag,n+1,n+1,"pair:noattraction_flag"); + memory->create(limit_damping,n+1,n+1,"pair:limit_damping"); onerad_dynamic = new double[n+1]; onerad_frozen = new double[n+1]; @@ -796,7 +796,7 @@ void PairGranular::coeff(int narg, char **arg) roll_model_one = ROLL_NONE; twist_model_one = TWIST_NONE; damping_model_one = VISCOELASTIC; - int na_flag = 0; + int ld_flag = 0; int iarg = 2; while (iarg < narg) { @@ -971,8 +971,8 @@ void PairGranular::coeff(int narg, char **arg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); cutoff_one = utils::numeric(FLERR,arg[iarg+1],false,lmp); iarg += 2; - } else if (strcmp(arg[iarg], "no_attraction") == 0) { - na_flag = 1; + } else if (strcmp(arg[iarg], "limit_damping") == 0) { + ld_flag = 1; iarg += 1; } else error->all(FLERR, "Illegal pair coeff command"); } @@ -991,11 +991,11 @@ void PairGranular::coeff(int narg, char **arg) 27.467*powint(cor,4)-18.022*powint(cor,5)+4.8218*powint(cor,6); } else damp = normal_coeffs_one[1]; - if(na_flag and normal_model_one == JKR) - error->all(FLERR,"Cannot turn off attraction with JKR pairstyle"); + if(ld_flag and normal_model_one == JKR) + error->all(FLERR,"Cannot limit damping with JKR model"); - if(na_flag and normal_model_one == DMT) - error->all(FLERR,"Cannot turn off attraction with DMT pairstyle"); + if(ld_flag and normal_model_one == DMT) + error->all(FLERR,"Cannot limit damping with DMT model"); for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { @@ -1039,7 +1039,7 @@ void PairGranular::coeff(int narg, char **arg) cutoff_type[i][j] = cutoff_type[j][i] = cutoff_one; - noattraction_flag[i][j] = na_flag; + limit_damping[i][j] = ld_flag; setflag[i][j] = 1; count++; @@ -1319,7 +1319,7 @@ void PairGranular::write_restart(FILE *fp) fwrite(&tangential_model[i][j],sizeof(int),1,fp); fwrite(&roll_model[i][j],sizeof(int),1,fp); fwrite(&twist_model[i][j],sizeof(int),1,fp); - fwrite(&noattraction_flag[i][j],sizeof(int),1,fp); + fwrite(&limit_damping[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); @@ -1350,7 +1350,7 @@ void PairGranular::read_restart(FILE *fp) utils::sfread(FLERR,&tangential_model[i][j],sizeof(int),1,fp,nullptr,error); utils::sfread(FLERR,&roll_model[i][j],sizeof(int),1,fp,nullptr,error); utils::sfread(FLERR,&twist_model[i][j],sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&noattraction_flag[i][j],sizeof(int),1,fp,nullptr,error); + utils::sfread(FLERR,&limit_damping[i][j],sizeof(int),1,fp,nullptr,error); utils::sfread(FLERR,normal_coeffs[i][j],sizeof(double),4,fp,nullptr,error); utils::sfread(FLERR,tangential_coeffs[i][j],sizeof(double),3,fp,nullptr,error); utils::sfread(FLERR,roll_coeffs[i][j],sizeof(double),3,fp,nullptr,error); @@ -1362,7 +1362,7 @@ void PairGranular::read_restart(FILE *fp) MPI_Bcast(&tangential_model[i][j],1,MPI_INT,0,world); MPI_Bcast(&roll_model[i][j],1,MPI_INT,0,world); MPI_Bcast(&twist_model[i][j],1,MPI_INT,0,world); - MPI_Bcast(&noattraction_flag[i][j],1,MPI_INT,0,world); + MPI_Bcast(&limit_damping[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); @@ -1548,7 +1548,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, Fdamp = -damp_normal_prefactor*vnnr; Fntot = Fne + Fdamp; - if(noattraction_flag[itype][jtype] and Fntot < 0.0) Fntot = 0.0; + if(limit_damping[itype][jtype] and Fntot < 0.0) Fntot = 0.0; jnum = list->numneigh[i]; jlist = list->firstneigh[i]; diff --git a/src/GRANULAR/pair_granular.h b/src/GRANULAR/pair_granular.h index 1e4e8c8fc4..49eb457f0d 100644 --- a/src/GRANULAR/pair_granular.h +++ b/src/GRANULAR/pair_granular.h @@ -70,7 +70,7 @@ class PairGranular : public Pair { // model choices int **normal_model, **damping_model; int **tangential_model, **roll_model, **twist_model; - int **noattractionflag; + int **limit_damping; // history flags int normal_history, tangential_history, roll_history, twist_history; diff --git a/src/KOKKOS/pair_gran_hooke_history_kokkos.cpp b/src/KOKKOS/pair_gran_hooke_history_kokkos.cpp index caa19ffdef..8d853b7dae 100644 --- a/src/KOKKOS/pair_gran_hooke_history_kokkos.cpp +++ b/src/KOKKOS/pair_gran_hooke_history_kokkos.cpp @@ -392,7 +392,7 @@ void PairGranHookeHistoryKokkos::operator()(TagPairGranHookeHistoryC F_FLOAT damp = meff*gamman*vnnr*rsqinv; F_FLOAT ccel = kn*(radsum-r)*rinv - damp; - if(noattractionflag & ccel < 0.0) ccel = 0.0; + if(limit_damping & ccel < 0.0) ccel = 0.0; // relative velocities From 8b37f3c04448bc3d3ed7a9f416cde67e8ad8e802 Mon Sep 17 00:00:00 2001 From: Joel Clemmer Date: Wed, 31 Mar 2021 17:25:41 -0600 Subject: [PATCH 05/15] Fixing incorrect arg for pair gran --- doc/src/fix_wall_gran.rst | 14 ++++---------- doc/src/fix_wall_gran_region.rst | 9 ++------- doc/src/pair_gran.rst | 6 +++--- doc/src/pair_granular.rst | 2 +- src/GRANULAR/pair_gran_hooke_history.cpp | 6 ++---- 5 files changed, 12 insertions(+), 25 deletions(-) diff --git a/doc/src/fix_wall_gran.rst b/doc/src/fix_wall_gran.rst index bfe06bf7a0..02d52ef1a6 100644 --- a/doc/src/fix_wall_gran.rst +++ b/doc/src/fix_wall_gran.rst @@ -29,6 +29,7 @@ Syntax gamma_t = damping coefficient for collisions in tangential direction (1/time units or 1/time-distance units - see discussion below) xmu = static yield criterion (unitless value between 0.0 and 1.0e4) dampflag = 0 or 1 if tangential damping force is excluded or included + optional keyword = *limit_damping*, limit damping to prevent attractive interaction .. parsed-literal:: @@ -45,7 +46,7 @@ Syntax radius = cylinder radius (distance units) * zero or more keyword/value pairs may be appended to args -* keyword = *wiggle* or *shear* or *contacts* or *no_attraction* +* keyword = *wiggle* or *shear* or *contacts* .. parsed-literal:: @@ -58,8 +59,6 @@ Syntax vshear = magnitude of shear velocity (velocity units) *contacts* value = none generate contact information for each particle - *no_attraction* value = none - turn off possibility of attractive interactions Examples @@ -97,7 +96,8 @@ Specifically, delta = radius - r = overlap of particle with wall, m_eff = mass of particle, and the effective radius of contact = RiRj/Ri+Rj is set to the radius of the particle. -The parameters *Kn*\ , *Kt*\ , *gamma_n*, *gamma_t*, *xmu* and *dampflag* +The parameters *Kn*\ , *Kt*\ , *gamma_n*, *gamma_t*, *xmu*, *dampflag*, +and the optional keyword *limit_damping* have the same meaning and units as those specified with the :doc:`pair_style gran/\* ` commands. This means a NULL can be used for either *Kt* or *gamma_t* as described on that page. If a @@ -177,12 +177,6 @@ the clockwise direction for *vshear* > 0 or counter-clockwise for *vshear* < 0. In this case, *vshear* is the tangential velocity of the wall at whatever *radius* has been defined. -If a particle is moving away from the wall while in contact, there -is a possibility that the particle could experience an effective attractive -force due to damping. If the *no_attraction* keyword is used, this fix -will zero out the normal component of the force if there is an effective -attractive force. This keyword cannot be used with the JKR or DMT models. - """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" This fix writes the shear friction state of atoms interacting with the diff --git a/doc/src/fix_wall_gran_region.rst b/doc/src/fix_wall_gran_region.rst index 00eb19ba1c..b7a9ff6bd7 100644 --- a/doc/src/fix_wall_gran_region.rst +++ b/doc/src/fix_wall_gran_region.rst @@ -183,7 +183,8 @@ radius - r = overlap of particle with wall, m_eff = mass of particle, and the effective radius of contact is just the radius of the particle. -The parameters *Kn*\ , *Kt*\ , *gamma_n*, *gamma_t*, *xmu* and *dampflag* +The parameters *Kn*\ , *Kt*\ , *gamma_n*, *gamma_t*, *xmu*, *dampflag*, +and the optional keyword *limit_damping* have the same meaning and units as those specified with the :doc:`pair_style gran/\* ` commands. This means a NULL can be used for either *Kt* or *gamma_t* as described on that page. If a @@ -201,12 +202,6 @@ values for the 6 wall/particle coefficients than for particle/particle interactions. E.g. if you wish to model the wall as a different material. -If a particle is moving away from the wall while in contact, there -is a possibility that the particle could experience an effective attractive -force due to damping. If the *no_attraction* keyword is used, this fix -will zero out the normal component of the force if there is an effective -attractive force. This keyword cannot be used with the JKR or DMT models. - Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" diff --git a/doc/src/pair_gran.rst b/doc/src/pair_gran.rst index 8145f7fb7a..b918e38da6 100644 --- a/doc/src/pair_gran.rst +++ b/doc/src/pair_gran.rst @@ -40,8 +40,8 @@ Syntax .. parsed-literal:: - *no_attraction* value = none - turn off possibility of attractive interactions + *limit_damping* value = none + limit damping to prevent attractive interaction .. note:: @@ -217,7 +217,7 @@ potential. If two particles are moving away from each other while in contact, there is a possibility that the particles could experience an effective attractive -force due to damping. If the *no_attraction* keyword is used, this fix +force due to damping. If the *limit_damping* keyword is used, this fix will zero out the normal component of the force if there is an effective attractive force. diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index 0076994e25..04c9d45afd 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -625,7 +625,7 @@ Finally, the twisting torque on each particle is given by: If two particles are moving away from each other while in contact, there is a possibility that the particles could experience an effective attractive -force due to damping. If the optional *no_attraction* keyword is used, this fix +force due to damping. If the optional *limit_damping* keyword is used, this fix will zero out the normal component of the force if there is an effective attractive force. This keyword cannot be used with the JKR or DMT models. diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp index d3c1cf7672..a86af2116e 100644 --- a/src/GRANULAR/pair_gran_hooke_history.cpp +++ b/src/GRANULAR/pair_gran_hooke_history.cpp @@ -373,10 +373,8 @@ void PairGranHookeHistory::settings(int narg, char **arg) limit_damping = 0; if(narg == 7) { - if(strcmp(arg[6], "limit_damping")) - limit_damping = 1; - else - error->all(FLERR,"Illegal pair_style command"); + if(strcmp(arg[6], "limit_damping") == 0) limit_damping = 1; + else error->all(FLERR,"Illegal pair_style command"); } if (kn < 0.0 || kt < 0.0 || gamman < 0.0 || gammat < 0.0 || From b1faf17eeb4ee978e25a034f4804afdf176c2f46 Mon Sep 17 00:00:00 2001 From: Dan Bolintineanu Date: Sat, 3 Apr 2021 00:48:00 -0600 Subject: [PATCH 06/15] Updates fix wall/gran to mirror recent updates to pair granular. Also some minor changes related to the limit_damping option --- doc/src/pair_gran.rst | 4 +- doc/src/pair_granular.rst | 2 +- src/GRANULAR/fix_wall_gran.cpp | 174 ++- src/GRANULAR/pair_gran_hooke_history.cpp | 6 +- src/GRANULAR/pair_granular.cpp | 31 +- src/GRANULAR/pair_granular_corrected.cpp | 1821 ++++++++++++++++++++++ 6 files changed, 1968 insertions(+), 70 deletions(-) create mode 100644 src/GRANULAR/pair_granular_corrected.cpp diff --git a/doc/src/pair_gran.rst b/doc/src/pair_gran.rst index b918e38da6..7d6189fda1 100644 --- a/doc/src/pair_gran.rst +++ b/doc/src/pair_gran.rst @@ -36,7 +36,7 @@ Syntax * xmu = static yield criterion (unitless value between 0.0 and 1.0e4) * dampflag = 0 or 1 if tangential damping force is excluded or included -* keyword = *no_attraction* +* keyword = *limit_damping* .. parsed-literal:: @@ -61,6 +61,8 @@ Examples pair_style gran/hooke/history 200000.0 NULL 50.0 NULL 0.5 1 pair_style gran/hooke 200000.0 70000.0 50.0 30.0 0.5 0 + pair_style gran/hooke 200000.0 70000.0 50.0 30.0 0.5 0 limit_damping + Description """"""""""" diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index 04c9d45afd..e0a33f31d0 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -24,7 +24,7 @@ Examples pair_coeff * * hooke 1000.0 50.0 tangential linear_history 500.0 1.0 0.4 damping mass_velocity pair_style granular - pair_coeff * * hertz 1000.0 50.0 tangential mindlin 1000.0 1.0 0.4 + pair_coeff * * hertz 1000.0 50.0 tangential mindlin 1000.0 1.0 0.4 limit_damping pair_style granular pair_coeff * * hertz/material 1e8 0.3 0.3 tangential mindlin_rescale NULL 1.0 0.4 damping tsuji diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index c9d00ccdd3..1d27e21162 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -53,7 +53,8 @@ enum{NONE,CONSTANT,EQUAL}; enum {NORMAL_HOOKE, NORMAL_HERTZ, HERTZ_MATERIAL, DMT, JKR}; enum {VELOCITY, MASS_VELOCITY, VISCOELASTIC, TSUJI}; enum {TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY, - TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE}; + TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE, + TANGENTIAL_MINDLIN_FORCE, TANGENTIAL_MINDLIN_RESCALE_FORCE}; enum {TWIST_NONE, TWIST_SDS, TWIST_MARSHALL}; enum {ROLL_NONE, ROLL_SDS}; @@ -216,7 +217,9 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : iarg += 4; } else if ((strcmp(arg[iarg+1], "linear_history") == 0) || (strcmp(arg[iarg+1], "mindlin") == 0) || - (strcmp(arg[iarg+1], "mindlin_rescale") == 0)) { + (strcmp(arg[iarg+1], "mindlin_rescale") == 0) || + (strcmp(arg[iarg+1], "mindlin/force") == 0) || + (strcmp(arg[iarg+1], "mindlin_rescale/force") == 0)) { if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, " "not enough parameters provided for tangential model"); @@ -226,8 +229,14 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : tangential_model = TANGENTIAL_MINDLIN; else if (strcmp(arg[iarg+1], "mindlin_rescale") == 0) tangential_model = TANGENTIAL_MINDLIN_RESCALE; + else if (strcmp(arg[iarg+1], "mindlin/force") == 0) + tangential_model = TANGENTIAL_MINDLIN_FORCE; + else if (strcmp(arg[iarg+1], "mindlin_rescale/force") == 0) + tangential_model = TANGENTIAL_MINDLIN_RESCALE_FORCE; if ((tangential_model == TANGENTIAL_MINDLIN || - tangential_model == TANGENTIAL_MINDLIN_RESCALE) && + tangential_model == TANGENTIAL_MINDLIN_RESCALE || + tangential_model == TANGENTIAL_MINDLIN_FORCE || + tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) && (strcmp(arg[iarg+2], "NULL") == 0)) { if (normal_model == NORMAL_HERTZ || normal_model == NORMAL_HOOKE) { error->all(FLERR, "NULL setting for Mindlin tangential " @@ -306,14 +315,20 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : } } size_history = 3*tangential_history + 3*roll_history + twist_history; + //Unlike the pair style, the wall style does not have a 'touch' + //array. Hence, an additional entry in the history is used to + //determine if particles previously contacted for JKR cohesion purposes. if (normal_model == JKR) size_history += 1; - if (tangential_model == TANGENTIAL_MINDLIN_RESCALE) size_history += 1; + if (tangential_model == TANGENTIAL_MINDLIN_RESCALE || + tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) size_history += 1; } - if(normal_model == JKR) - error->all(FLERR,"Cannot limit damping with JRK model"); - if(normal_model == DMT) - error->all(FLERR,"Cannot limit damping with DMT model"); + if (limit_damping and normal_model == JKR) + error->all(FLERR,"Illegal pair_coeff command, " + "cannot limit damping with JRK model"); + if (limit_damping and normal_model == DMT) + error->all(FLERR,"Illegal pair_coeff command, " + "Cannot limit damping with DMT model"); // wallstyle args @@ -509,7 +524,8 @@ void FixWallGran::init() roll_history_index += 1; twist_history_index += 1; } - if (tangential_model == TANGENTIAL_MINDLIN_RESCALE) { + if (tangential_model == TANGENTIAL_MINDLIN_RESCALE || + tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) { roll_history_index += 1; twist_history_index += 1; } @@ -1120,7 +1136,8 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz, double signtwist, magtwist, magtortwist, Mtcrit; double tortwist1, tortwist2, tortwist3; - double shrmag,rsht; + double shrmag,rsht,prjmag; + bool frameupdate; r = sqrt(rsq); E = normal_coeffs[0]; @@ -1195,12 +1212,18 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz, Fdamp = -damp_normal_prefactor*vnnr; Fntot = Fne + Fdamp; - if(limit_damping and Fntot < 0.0) Fntot = 0.0; + if (limit_damping and Fntot < 0.0) Fntot = 0.0; //**************************************** // tangential force, including history effects //**************************************** + // For linear, mindlin, mindlin_rescale: + // history = cumulative tangential displacement + // + // For mindlin/force, mindlin_rescale/force: + // history = cumulative tangential elastic force + // tangential component vt1 = vr1 - vn1; vt2 = vr2 - vn2; @@ -1242,69 +1265,115 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz, int thist2 = thist1 + 1; if (tangential_history) { - if (tangential_model == TANGENTIAL_MINDLIN) { + if (tangential_model == TANGENTIAL_MINDLIN || + tangential_model == TANGENTIAL_MINDLIN_FORCE) { k_tangential *= a; } - else if (tangential_model == TANGENTIAL_MINDLIN_RESCALE) { + else if (tangential_model == + TANGENTIAL_MINDLIN_RESCALE || + tangential_model == + TANGENTIAL_MINDLIN_RESCALE_FORCE){ k_tangential *= a; - if (a < history[3]) { //On unloading, rescale the shear displacements + // on unloading, rescale the shear displacements/force + if (a < history[thist2+1]) { double factor = a/history[thist2+1]; history[thist0] *= factor; history[thist1] *= factor; history[thist2] *= factor; } } - shrmag = sqrt(history[thist0]*history[thist0] + - history[thist1]*history[thist1] + - history[thist2]*history[thist2]); + // rotate and update displacements. // see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 if (history_update) { rsht = history[thist0]*nx + history[thist1]*ny + history[thist2]*nz; - if (fabs(rsht) < EPSILON) rsht = 0; - if (rsht > 0) { - // if rhst == shrmag, contacting pair has rotated 90 deg in one step, - // in which case you deserve a crash! - scalefac = shrmag/(shrmag - rsht); + if (tangential_model == TANGENTIAL_MINDLIN_FORCE || + tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) + frameupdate = fabs(rsht) > EPSILON*Fscrit; + else + frameupdate = fabs(rsht)*k_tangential > EPSILON*Fscrit; + if (frameupdate) { + shrmag = sqrt(history[thist0]*history[thist0] + + history[thist1]*history[thist1] + + history[thist2]*history[thist2]); + // projection history[thist0] -= rsht*nx; history[thist1] -= rsht*ny; history[thist2] -= rsht*nz; + // also rescale to preserve magnitude + prjmag = sqrt(history[0]*history[0] + history[1]*history[1] + + history[2]*history[2]); + if (prjmag > 0) scalefac = shrmag/prjmag; + else scalefac = 0; history[thist0] *= scalefac; history[thist1] *= scalefac; history[thist2] *= scalefac; } // update history + if (tangential_model == TANGENTIAL_HISTORY || + tangential_model == TANGENTIAL_MINDLIN || + tangential_model == TANGENTIAL_MINDLIN_RESCALE) { history[thist0] += vtr1*dt; history[thist1] += vtr2*dt; history[thist2] += vtr3*dt; + } else{ + // tangential force + // see e.g. eq. 18 of Thornton et al, Pow. Tech. 2013, v223,p30-46 + history[thist0] -= k_tangential*vtr1*dt; + history[thist1] -= k_tangential*vtr2*dt; + history[thist2] -= k_tangential*vtr3*dt; + } + if (tangential_model == TANGENTIAL_MINDLIN_RESCALE || + tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) + history[thist2+1] = a; } // tangential forces = history + tangential velocity damping - fs1 = -k_tangential*history[thist0] - damp_tangential*vtr1; - fs2 = -k_tangential*history[thist1] - damp_tangential*vtr2; - fs3 = -k_tangential*history[thist2] - damp_tangential*vtr3; + if (tangential_model == TANGENTIAL_HISTORY || + tangential_model == TANGENTIAL_MINDLIN || + tangential_model == TANGENTIAL_MINDLIN_RESCALE) { + fs1 = -k_tangential*history[thist0] - damp_tangential*vtr1; + fs2 = -k_tangential*history[thist1] - damp_tangential*vtr2; + fs3 = -k_tangential*history[thist2] - damp_tangential*vtr3; + } else { + fs1 = history[thist0] - damp_tangential*vtr1; + fs2 = history[thist1] - damp_tangential*vtr2; + fs3 = history[thist2] - damp_tangential*vtr3; + } // rescale frictional displacements and forces if needed Fscrit = tangential_coeffs[2] * Fncrit; fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); if (fs > Fscrit) { + shrmag = sqrt(history[thist0]*history[thist0] + + history[thist1]*history[thist1] + + history[thist2]*history[thist2]); if (shrmag != 0.0) { - history[thist0] = -1.0/k_tangential*(Fscrit*fs1/fs + - damp_tangential*vtr1); - history[thist1] = -1.0/k_tangential*(Fscrit*fs2/fs + - damp_tangential*vtr2); - history[thist2] = -1.0/k_tangential*(Fscrit*fs3/fs + - damp_tangential*vtr3); + if (tangential_model == TANGENTIAL_HISTORY || + tangential_model == TANGENTIAL_MINDLIN || + tangential_model == + TANGENTIAL_MINDLIN_RESCALE) { + history[thist0] = -1.0/k_tangential*(Fscrit*fs1/fs + + damp_tangential*vtr1); + history[thist1] = -1.0/k_tangential*(Fscrit*fs2/fs + + damp_tangential*vtr2); + history[thist2] = -1.0/k_tangential*(Fscrit*fs3/fs + + damp_tangential*vtr3); + } else { + history[thist0] = Fscrit*fs1/fs + damp_tangential*vtr1; + history[thist1] = Fscrit*fs2/fs + damp_tangential*vtr2; + history[thist2] = 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; + fs = damp_tangential*vrel; + if (vrel != 0.0) Ft = MIN(Fscrit,fs) / vrel; else Ft = 0.0; fs1 = -Ft*vtr1; fs2 = -Ft*vtr2; @@ -1315,14 +1384,14 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz, // rolling resistance //**************************************** - if (roll_model != ROLL_NONE || twist_model != NONE) { + if (roll_model != ROLL_NONE || twist_model != TWIST_NONE) { relrot1 = omega[0]; relrot2 = omega[1]; relrot3 = omega[2]; } if (roll_model != ROLL_NONE) { - - // rolling velocity, see eq. 31 of Wang et al, Particuology v 23, p 49 (2015) + // 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) @@ -1334,21 +1403,29 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz, 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[0]; + damp_roll = roll_coeffs[1]; + Frcrit = roll_coeffs[2] * Fncrit; if (history_update) { - if (fabs(rolldotn) < EPSILON) rolldotn = 0; - if (rolldotn > 0) { // rotate into tangential plane - scalefac = rollmag/(rollmag - rolldotn); + rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz; + frameupdate = fabs(rolldotn)*k_roll > EPSILON*Frcrit; + if (frameupdate) { // rotate into tangential plane + rollmag = sqrt(history[rhist0]*history[rhist0] + + history[rhist1]*history[rhist1] + + history[rhist2]*history[rhist2]); + // projection history[rhist0] -= rolldotn*nx; history[rhist1] -= rolldotn*ny; history[rhist2] -= rolldotn*nz; + // also rescale to preserve magnitude + prjmag = sqrt(history[rhist0]*history[rhist0] + + history[rhist1]*history[rhist1] + + history[rhist2]*history[rhist2]); + + if (prjmag > 0) scalefac = rollmag/prjmag; + else scalefac = 0; history[rhist0] *= scalefac; history[rhist1] *= scalefac; history[rhist2] *= scalefac; @@ -1358,17 +1435,16 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz, history[rhist2] += vrl3*dt; } - k_roll = roll_coeffs[0]; - damp_roll = roll_coeffs[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[2] * Fncrit; - fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3); if (fr > Frcrit) { + rollmag = sqrt(history[rhist0]*history[rhist0] + + history[rhist1]*history[rhist1] + + history[rhist2]*history[rhist2]); 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); diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp index a86af2116e..4f7a6e20ec 100644 --- a/src/GRANULAR/pair_gran_hooke_history.cpp +++ b/src/GRANULAR/pair_gran_hooke_history.cpp @@ -237,7 +237,7 @@ void PairGranHookeHistory::compute(int eflag, int vflag) damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; - if(limit_damping and ccel < 0.0) ccel = 0.0; + if (limit_damping and ccel < 0.0) ccel = 0.0; // relative velocities @@ -372,8 +372,8 @@ void PairGranHookeHistory::settings(int narg, char **arg) if (dampflag == 0) gammat = 0.0; limit_damping = 0; - if(narg == 7) { - if(strcmp(arg[6], "limit_damping") == 0) limit_damping = 1; + if (narg == 7) { + if (strcmp(arg[6], "limit_damping") == 0) limit_damping = 1; else error->all(FLERR,"Illegal pair_style command"); } diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index dc829879e8..8bf9b0c4ed 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -365,13 +365,13 @@ void PairGranular::compute(int eflag, int vflag) damp_normal = a*meff; } else if (damping_model[itype][jtype] == TSUJI) { damp_normal = sqrt(meff*knfac); - } + } else damp_normal = 0.0; damp_normal_prefactor = normal_coeffs[itype][jtype][1]*damp_normal; Fdamp = -damp_normal_prefactor*vnnr; Fntot = Fne + Fdamp; - if(limit_damping[itype][jtype] and Fntot < 0.0) Fntot = 0.0; + if (limit_damping[itype][jtype] and Fntot < 0.0) Fntot = 0.0; //**************************************** // tangential force, including history effects @@ -542,20 +542,17 @@ void PairGranular::compute(int eflag, int vflag) 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) - // - 0.5*((radj-radi)/radsum)*vtr1; - // - 0.5*((radj-radi)/radsum)*vtr2; - // - 0.5*((radj-radi)/radsum)*vtr3; } //**************************************** // rolling resistance //**************************************** if (roll_model[itype][jtype] != ROLL_NONE) { + // 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); vrl2 = Reff*(relrot3*nx - relrot1*nz); vrl3 = Reff*(relrot1*ny - relrot2*nx); @@ -974,12 +971,12 @@ void PairGranular::coeff(int narg, char **arg) } else if (strcmp(arg[iarg], "limit_damping") == 0) { ld_flag = 1; iarg += 1; - } else error->all(FLERR, "Illegal pair coeff command"); + } else error->all(FLERR, "Illegal pair_coeff command"); } // error not to specify normal or tangential model if ((normal_model_one < 0) || (tangential_model_one < 0)) - error->all(FLERR, "Illegal pair coeff command, " + error->all(FLERR, "Illegal pair_coeff command, " "must specify normal or tangential contact model"); int count = 0; @@ -991,11 +988,13 @@ void PairGranular::coeff(int narg, char **arg) 27.467*powint(cor,4)-18.022*powint(cor,5)+4.8218*powint(cor,6); } else damp = normal_coeffs_one[1]; - if(ld_flag and normal_model_one == JKR) - error->all(FLERR,"Cannot limit damping with JKR model"); + if (ld_flag and normal_model_one == JKR) + error->all(FLERR,"Illegal pair_coeff command, " + "Cannot limit damping with JKR model"); - if(ld_flag and normal_model_one == DMT) - error->all(FLERR,"Cannot limit damping with DMT model"); + if (ld_flag and normal_model_one == DMT) + error->all(FLERR,"Illegal pair_coeff command, " + "Cannot limit damping with DMT model"); for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { diff --git a/src/GRANULAR/pair_granular_corrected.cpp b/src/GRANULAR/pair_granular_corrected.cpp new file mode 100644 index 0000000000..bcf262aa2c --- /dev/null +++ b/src/GRANULAR/pair_granular_corrected.cpp @@ -0,0 +1,1821 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://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 "pair_granular.h" + +#include +#include + +#include "atom.h" +#include "force.h" +#include "update.h" +#include "modify.h" +#include "fix.h" +#include "fix_dummy.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 "math_special.h" + + +using namespace LAMMPS_NS; +using namespace MathConst; +using namespace MathSpecial; + +#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 4.0/3.0 // 4/3 +#define THREEQUARTERS 0.75 // 3/4 + +#define EPSILON 1e-10 + +enum {HOOKE, HERTZ, HERTZ_MATERIAL, DMT, JKR}; +enum {VELOCITY, MASS_VELOCITY, VISCOELASTIC, TSUJI}; +enum {TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY, + TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE, + TANGENTIAL_MINDLIN_FORCE, TANGENTIAL_MINDLIN_RESCALE_FORCE}; +enum {TWIST_NONE, TWIST_SDS, TWIST_MARSHALL}; +enum {ROLL_NONE, ROLL_SDS}; + +/* ---------------------------------------------------------------------- */ + +PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp) +{ + single_enable = 1; + no_virial_fdotr_compute = 1; + centroidstressflag = CENTROID_NOTAVAIL; + + single_extra = 12; + svector = new double[single_extra]; + + neighprev = 0; + + nmax = 0; + mass_rigid = nullptr; + + onerad_dynamic = nullptr; + onerad_frozen = nullptr; + maxrad_dynamic = nullptr; + maxrad_frozen = nullptr; + + history_transfer_factors = nullptr; + + 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; + + // create dummy fix as placeholder for FixNeighHistory + // this is so final order of Modify:fix will conform to input script + + fix_history = nullptr; + modify->add_fix("NEIGH_HISTORY_GRANULAR_DUMMY all DUMMY"); + fix_dummy = (FixDummy *) modify->fix[modify->nfix-1]; +} + +/* ---------------------------------------------------------------------- */ + +PairGranular::~PairGranular() +{ + delete [] svector; + + if (!fix_history) modify->delete_fix("NEIGH_HISTORY_GRANULAR_DUMMY"); + else modify->delete_fix("NEIGH_HISTORY_GRANULAR"); + + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + memory->destroy(cutoff_type); + + memory->destroy(normal_coeffs); + memory->destroy(tangential_coeffs); + memory->destroy(roll_coeffs); + memory->destroy(twist_coeffs); + + memory->destroy(Emod); + memory->destroy(poiss); + + memory->destroy(normal_model); + memory->destroy(damping_model); + memory->destroy(tangential_model); + memory->destroy(roll_model); + memory->destroy(twist_model); + + delete [] onerad_dynamic; + delete [] onerad_frozen; + delete [] maxrad_dynamic; + delete [] maxrad_frozen; + } + + memory->destroy(mass_rigid); +} + +/* ---------------------------------------------------------------------- */ + +void PairGranular::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; + double Reff, delta, dR, dR2, dist_to_contact; + + double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; + double wr1,wr2,wr3; + double vtr1,vtr2,vtr3,vrel; + + double knfac, damp_normal=0.0, damp_normal_prefactor; + double k_tangential, damp_tangential; + double Fne, Ft, Fdamp, Fntot, Fncrit, Fscrit, Frcrit; + double fs, fs1, fs2, fs3, tor1, tor2, tor3; + + double mi,mj,meff; + double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3; + + // 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; + + // rolling + double k_roll, damp_roll; + double 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,prjmag; + bool frameupdate; + int *ilist,*jlist,*numneigh,**firstneigh; + int *touch,**firsttouch; + double *history,*allhistory,**firsthistory; + + bool touchflag = false; + const bool historyupdate = (update->setupflag) ? false : true; + + ev_init(eflag,vflag); + + // 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; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + if (use_history) { + 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]; + if (use_history) { + 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/radsum; + touchflag = false; + + if (normal_model[itype][jtype] == JKR) { + E *= THREEQUARTERS; + if (touch[jj]) { + R2 = Reff*Reff; + coh = normal_coeffs[itype][jtype][3]; + a = cbrt(9.0*MY_PI*coh*R2/(4*E)); + delta_pulloff = a*a/Reff - 2*sqrt(MY_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 + if (use_history) { + 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_model[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; + // in case sqrt(0) < 0 due to precision issues + sqrt1 = MAX(0, t0*(t1+2*t2)); + t4 = cbrt(t1+t2+THREEROOT3*MY_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*MY_PI*R2/(E*t6)); + a = INVROOT6*(t6 + sqrt(sqrt3)); + a2 = a*a; + knfac = normal_coeffs[itype][jtype][0]*a; + Fne = knfac*a2/Reff - MY_2PI*a2*sqrt(4*coh*E/(MY_PI*a)); + } else { + knfac = E; // Hooke + Fne = knfac*delta; + a = sqrt(dR); + if (normal_model[itype][jtype] != HOOKE) { + Fne *= a; + knfac *= a; + } + if (normal_model[itype][jtype] == DMT) + Fne -= 4*MY_PI*normal_coeffs[itype][jtype][3]*Reff; + } + + // NOTE: consider restricting Hooke to only have + // 'velocity' as an option for damping? + + if (damping_model[itype][jtype] == VELOCITY) { + damp_normal = 1; + } else if (damping_model[itype][jtype] == MASS_VELOCITY) { + damp_normal = meff; + } else if (damping_model[itype][jtype] == VISCOELASTIC) { + damp_normal = a*meff; + } else if (damping_model[itype][jtype] == TSUJI) { + damp_normal = sqrt(meff*knfac); + } + + damp_normal_prefactor = normal_coeffs[itype][jtype][1]*damp_normal; + Fdamp = -damp_normal_prefactor*vnnr; + + Fntot = Fne + Fdamp; + + //**************************************** + // tangential force, including history effects + //**************************************** + + // For linear, mindlin, mindlin_rescale: + // history = cumulative tangential displacement + // + // For mindlin/force, mindlin_rescale/force: + // history = cumulative tangential elastic force + + // 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_model[itype][jtype] == JKR) { + F_pulloff = 3*MY_PI*coh*Reff; + Fncrit = fabs(Fne + 2*F_pulloff); + } else if (normal_model[itype][jtype] == DMT) { + F_pulloff = 4*MY_PI*coh*Reff; + Fncrit = fabs(Fne + 2*F_pulloff); + } else { + Fncrit = fabs(Fntot); + } + Fscrit = tangential_coeffs[itype][jtype][2] * Fncrit; + + //------------------------------ + // tangential forces + //------------------------------ + k_tangential = tangential_coeffs[itype][jtype][0]; + damp_tangential = tangential_coeffs[itype][jtype][1] * + damp_normal_prefactor; + + if (tangential_history) { + if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN || + tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_FORCE) { + k_tangential *= a; + } else if (tangential_model[itype][jtype] == + TANGENTIAL_MINDLIN_RESCALE || + tangential_model[itype][jtype] == + TANGENTIAL_MINDLIN_RESCALE_FORCE) { + k_tangential *= a; + // on unloading, rescale the shear displacements/force + if (a < history[3]) { + double factor = a/history[3]; + history[0] *= factor; + history[1] *= factor; + history[2] *= factor; + } + } + // rotate and update displacements / force. + // 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 (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_FORCE || + tangential_model[itype][jtype] == + TANGENTIAL_MINDLIN_RESCALE_FORCE) + frameupdate = fabs(rsht) > EPSILON*Fscrit; + else + frameupdate = fabs(rsht)*k_tangential > EPSILON*Fscrit; + if (frameupdate) { + shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + + history[2]*history[2]); + // projection + history[0] -= rsht*nx; + history[1] -= rsht*ny; + history[2] -= rsht*nz; + + // also rescale to preserve magnitude + prjmag = sqrt(history[0]*history[0] + history[1]*history[1] + + history[2]*history[2]); + if (prjmag > 0) scalefac = shrmag/prjmag; + else scalefac = 0; + history[0] *= scalefac; + history[1] *= scalefac; + history[2] *= scalefac; + } + // update history + if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY || + tangential_model[itype][jtype] == TANGENTIAL_MINDLIN || + tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) { + // tangential displacement + history[0] += vtr1*dt; + history[1] += vtr2*dt; + history[2] += vtr3*dt; + } else { + // tangential force + // see e.g. eq. 18 of Thornton et al, Pow. Tech. 2013, v223,p30-46 + history[0] -= k_tangential*vtr1*dt; + history[1] -= k_tangential*vtr2*dt; + history[2] -= k_tangential*vtr3*dt; + } + if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE || + tangential_model[itype][jtype] == + TANGENTIAL_MINDLIN_RESCALE_FORCE) + history[3] = a; + } + + // tangential forces = history + tangential velocity damping + if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY || + tangential_model[itype][jtype] == TANGENTIAL_MINDLIN || + tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) { + fs1 = -k_tangential*history[0] - damp_tangential*vtr1; + fs2 = -k_tangential*history[1] - damp_tangential*vtr2; + fs3 = -k_tangential*history[2] - damp_tangential*vtr3; + } else { + fs1 = history[0] - damp_tangential*vtr1; + fs2 = history[1] - damp_tangential*vtr2; + fs3 = history[2] - damp_tangential*vtr3; + } + + // rescale frictional displacements and forces if needed + fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); + if (fs > Fscrit) { + shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + + history[2]*history[2]); + if (shrmag != 0.0) { + if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY || + tangential_model[itype][jtype] == TANGENTIAL_MINDLIN || + tangential_model[itype][jtype] == + TANGENTIAL_MINDLIN_RESCALE) { + 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); + } else { + history[0] = Fscrit*fs1/fs + damp_tangential*vtr1; + history[1] = Fscrit*fs2/fs + damp_tangential*vtr2; + history[2] = 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 = damp_tangential*vrel; + if (vrel != 0.0) Ft = MIN(Fscrit,fs) / vrel; + else Ft = 0.0; + fs1 = -Ft*vtr1; + fs2 = -Ft*vtr2; + fs3 = -Ft*vtr3; + } + + if (roll_model[itype][jtype] != ROLL_NONE || + twist_model[itype][jtype] != TWIST_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) + // - 0.5*((radj-radi)/radsum)*vtr1; + // - 0.5*((radj-radi)/radsum)*vtr2; + // - 0.5*((radj-radi)/radsum)*vtr3; + } + //**************************************** + // rolling resistance + //**************************************** + + if (roll_model[itype][jtype] != ROLL_NONE) { + vrl1 = Reff*(relrot2*nz - relrot3*ny); + vrl2 = Reff*(relrot3*nx - relrot1*nz); + vrl3 = Reff*(relrot1*ny - relrot2*nx); + + int rhist0 = roll_history_index; + int rhist1 = rhist0 + 1; + int rhist2 = rhist1 + 1; + + k_roll = roll_coeffs[itype][jtype][0]; + damp_roll = roll_coeffs[itype][jtype][1]; + Frcrit = roll_coeffs[itype][jtype][2] * Fncrit; + + if (historyupdate) { + rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz; + frameupdate = fabs(rolldotn)*k_roll > EPSILON*Frcrit; + if (frameupdate) { // rotate into tangential plane + rollmag = sqrt(history[rhist0]*history[rhist0] + + history[rhist1]*history[rhist1] + + history[rhist2]*history[rhist2]); + // projection + history[rhist0] -= rolldotn*nx; + history[rhist1] -= rolldotn*ny; + history[rhist2] -= rolldotn*nz; + // also rescale to preserve magnitude + prjmag = sqrt(history[rhist0]*history[rhist0] + + history[rhist1]*history[rhist1] + + history[rhist2]*history[rhist2]); + if (prjmag > 0) scalefac = rollmag/prjmag; + else scalefac = 0; + history[rhist0] *= scalefac; + history[rhist1] *= scalefac; + history[rhist2] *= scalefac; + } + history[rhist0] += vrl1*dt; + history[rhist1] += vrl2*dt; + history[rhist2] += vrl3*dt; + } + + 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 + + fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3); + if (fr > Frcrit) { + rollmag = sqrt(history[rhist0]*history[rhist0] + + history[rhist1]*history[rhist1] + + history[rhist2]*history[rhist2]); + 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; + } + } + + //**************************************** + // twisting torque, including history effects + //**************************************** + + if (twist_model[itype][jtype] != TWIST_NONE) { + // omega_T (eq 29 of Marshall) + magtwist = relrot1*nx + relrot2*ny + relrot3*nz; + if (twist_model[itype][jtype] == TWIST_MARSHALL) { + k_twist = 0.5*k_tangential*a*a;; // eq 32 of Marshall paper + damp_twist = 0.5*damp_tangential*a*a; + mu_twist = TWOTHIRDS*a*tangential_coeffs[itype][jtype][2]; + } else { + k_twist = twist_coeffs[itype][jtype][0]; + damp_twist = twist_coeffs[itype][jtype][1]; + mu_twist = twist_coeffs[itype][jtype][2]; + } + 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 = mu_twist*Fncrit; // 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 + } + } + + // 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; + + dist_to_contact = radi-0.5*delta; + torque[i][0] -= dist_to_contact*tor1; + torque[i][1] -= dist_to_contact*tor2; + torque[i][2] -= dist_to_contact*tor3; + + if (twist_model[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_model[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; + + dist_to_contact = radj-0.5*delta; + torque[j][0] -= dist_to_contact*tor1; + torque[j][1] -= dist_to_contact*tor2; + torque[j][2] -= dist_to_contact*tor3; + + if (twist_model[itype][jtype] != TWIST_NONE) { + torque[j][0] -= tortwist1; + torque[j][1] -= tortwist2; + torque[j][2] -= tortwist3; + } + if (roll_model[itype][jtype] != ROLL_NONE) { + torque[j][0] -= torroll1; + torque[j][1] -= torroll2; + torque[j][2] -= torroll3; + } + } + if (evflag) ev_tally_xyz(i,j,nlocal,force->newton_pair, + 0.0,0.0,fx,fy,fz,delx,dely,delz); + } + } + } +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairGranular::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(cutoff_type,n+1,n+1,"pair:cutoff_type"); + 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(Emod,n+1,n+1,"pair:Emod"); + memory->create(poiss,n+1,n+1,"pair:poiss"); + + memory->create(normal_model,n+1,n+1,"pair:normal_model"); + memory->create(damping_model,n+1,n+1,"pair:damping_model"); + memory->create(tangential_model,n+1,n+1,"pair:tangential_model"); + memory->create(roll_model,n+1,n+1,"pair:roll_model"); + memory->create(twist_model,n+1,n+1,"pair:twist_model"); + + 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 PairGranular::settings(int narg, char **arg) +{ + if (narg == 1) { + cutoff_global = utils::numeric(FLERR,arg[0],false,lmp); + } else { + cutoff_global = -1; // will be set based on particle sizes, model choice + } + + normal_history = tangential_history = 0; + roll_history = twist_history = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairGranular::coeff(int narg, char **arg) +{ + int normal_model_one, damping_model_one; + int tangential_model_one, roll_model_one, twist_model_one; + + double normal_coeffs_one[4]; + double tangential_coeffs_one[4]; + double roll_coeffs_one[4]; + double twist_coeffs_one[4]; + + double cutoff_one = -1; + + if (narg < 2) + error->all(FLERR,"Incorrect args for pair coefficients"); + + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); + utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); + + //Defaults + normal_model_one = tangential_model_one = -1; + roll_model_one = ROLL_NONE; + twist_model_one = TWIST_NONE; + damping_model_one = 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_model_one = HOOKE; + normal_coeffs_one[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp); // kn + normal_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp); // damping + iarg += 3; + } else if (strcmp(arg[iarg], "hertz") == 0) { + if (iarg + 2 >= narg) + error->all(FLERR,"Illegal pair_coeff command, " + "not enough parameters provided for Hertz option"); + normal_model_one = HERTZ; + normal_coeffs_one[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp); // kn + normal_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp); // damping + iarg += 3; + } else if (strcmp(arg[iarg], "hertz/material") == 0) { + if (iarg + 3 >= narg) + error->all(FLERR,"Illegal pair_coeff command, " + "not enough parameters provided for Hertz/material option"); + normal_model_one = HERTZ_MATERIAL; + normal_coeffs_one[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp); // E + normal_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp); // damping + normal_coeffs_one[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp); // Poisson's ratio + iarg += 4; + } 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_model_one = DMT; + normal_coeffs_one[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp); // E + normal_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp); // damping + normal_coeffs_one[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp); // Poisson's ratio + normal_coeffs_one[3] = utils::numeric(FLERR,arg[iarg+4],false,lmp); // 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_model_one = JKR; + normal_coeffs_one[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp); // E + normal_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp); // damping + normal_coeffs_one[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp); // Poisson's ratio + normal_coeffs_one[3] = utils::numeric(FLERR,arg[iarg+4],false,lmp); // cohesion + iarg += 5; + } else if (strcmp(arg[iarg], "damping") == 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_model_one = VELOCITY; + iarg += 1; + } else if (strcmp(arg[iarg+1], "mass_velocity") == 0) { + damping_model_one = MASS_VELOCITY; + iarg += 1; + } else if (strcmp(arg[iarg+1], "viscoelastic") == 0) { + damping_model_one = VISCOELASTIC; + iarg += 1; + } else if (strcmp(arg[iarg+1], "tsuji") == 0) { + damping_model_one = TSUJI; + iarg += 1; + } else error->all(FLERR, "Illegal pair_coeff command, " + "unrecognized damping model"); + iarg += 1; + } else if (strcmp(arg[iarg], "tangential") == 0) { + if (iarg + 1 >= narg) + error->all(FLERR,"Illegal pair_coeff command, must specify " + "tangential model after tangential keyword"); + if (strcmp(arg[iarg+1], "linear_nohistory") == 0) { + if (iarg + 3 >= narg) + error->all(FLERR,"Illegal pair_coeff command, " + "not enough parameters provided for tangential model"); + tangential_model_one = TANGENTIAL_NOHISTORY; + tangential_coeffs_one[0] = 0; + // gammat and friction coeff + tangential_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp); + tangential_coeffs_one[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp); + iarg += 4; + } else if ((strcmp(arg[iarg+1], "linear_history") == 0) || + (strcmp(arg[iarg+1], "mindlin") == 0) || + (strcmp(arg[iarg+1], "mindlin_rescale") == 0) || + (strcmp(arg[iarg+1], "mindlin/force") == 0) || + (strcmp(arg[iarg+1], "mindlin_rescale/force") == 0)) { + if (iarg + 4 >= narg) + error->all(FLERR,"Illegal pair_coeff command, " + "not enough parameters provided for tangential model"); + if (strcmp(arg[iarg+1], "linear_history") == 0) + tangential_model_one = TANGENTIAL_HISTORY; + else if (strcmp(arg[iarg+1], "mindlin") == 0) + tangential_model_one = TANGENTIAL_MINDLIN; + else if (strcmp(arg[iarg+1], "mindlin_rescale") == 0) + tangential_model_one = TANGENTIAL_MINDLIN_RESCALE; + else if (strcmp(arg[iarg+1], "mindlin/force") == 0) + tangential_model_one = TANGENTIAL_MINDLIN_FORCE; + else if (strcmp(arg[iarg+1], "mindlin_rescale/force") == 0) + tangential_model_one = TANGENTIAL_MINDLIN_RESCALE_FORCE; + tangential_history = 1; + if ((tangential_model_one == TANGENTIAL_MINDLIN || + tangential_model_one == TANGENTIAL_MINDLIN_RESCALE || + tangential_model_one == TANGENTIAL_MINDLIN_FORCE || + tangential_model_one == TANGENTIAL_MINDLIN_RESCALE_FORCE) && + (strcmp(arg[iarg+2], "NULL") == 0)) { + if (normal_model_one == HERTZ || normal_model_one == HOOKE) { + error->all(FLERR, "NULL setting for Mindlin tangential " + "stiffness requires a normal contact model that " + "specifies material properties"); + } + tangential_coeffs_one[0] = -1; + } else { + tangential_coeffs_one[0] = utils::numeric(FLERR,arg[iarg+2],false,lmp); // kt + } + // gammat and friction coeff + tangential_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+3],false,lmp); + tangential_coeffs_one[2] = utils::numeric(FLERR,arg[iarg+4],false,lmp); + iarg += 5; + } else { + error->all(FLERR, "Illegal pair_coeff command, " + "tangential model not recognized"); + } + } 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_model_one = ROLL_NONE; + iarg += 2; + } else if (strcmp(arg[iarg+1], "sds") == 0) { + if (iarg + 4 >= narg) + error->all(FLERR,"Illegal pair_coeff command, " + "not enough parameters provided for rolling model"); + roll_model_one = ROLL_SDS; + roll_history = 1; + // kR and gammaR and rolling friction coeff + roll_coeffs_one[0] = utils::numeric(FLERR,arg[iarg+2],false,lmp); + roll_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+3],false,lmp); + roll_coeffs_one[2] = utils::numeric(FLERR,arg[iarg+4],false,lmp); + iarg += 5; + } else { + error->all(FLERR, "Illegal pair_coeff command, " + "rolling friction model not recognized"); + } + } 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_model_one = TWIST_NONE; + iarg += 2; + } else if (strcmp(arg[iarg+1], "marshall") == 0) { + twist_model_one = TWIST_MARSHALL; + twist_history = 1; + iarg += 2; + } else if (strcmp(arg[iarg+1], "sds") == 0) { + if (iarg + 4 >= narg) + error->all(FLERR,"Illegal pair_coeff command, " + "not enough parameters provided for twist model"); + twist_model_one = TWIST_SDS; + twist_history = 1; + // kt and gammat and friction coeff + twist_coeffs_one[0] = utils::numeric(FLERR,arg[iarg+2],false,lmp); + twist_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+3],false,lmp); + twist_coeffs_one[2] = utils::numeric(FLERR,arg[iarg+4],false,lmp); + iarg += 5; + } else { + error->all(FLERR, "Illegal pair_coeff command, " + "twisting friction model not recognized"); + } + } else if (strcmp(arg[iarg], "cutoff") == 0) { + if (iarg + 1 >= narg) + error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); + cutoff_one = utils::numeric(FLERR,arg[iarg+1],false,lmp); + iarg += 2; + } else error->all(FLERR, "Illegal pair coeff command"); + } + + // error not to specify normal or tangential model + if ((normal_model_one < 0) || (tangential_model_one < 0)) + error->all(FLERR, "Illegal pair coeff command, " + "must specify normal or tangential contact model"); + + int count = 0; + double damp; + if (damping_model_one == TSUJI) { + double cor; + cor = normal_coeffs_one[1]; + damp = 1.2728-4.2783*cor+11.087*square(cor)-22.348*cube(cor)+ + 27.467*powint(cor,4)-18.022*powint(cor,5)+4.8218*powint(cor,6); + } else damp = normal_coeffs_one[1]; + + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + normal_model[i][j] = normal_model[j][i] = normal_model_one; + normal_coeffs[i][j][1] = normal_coeffs[j][i][1] = damp; + if (normal_model_one != HERTZ && normal_model_one != HOOKE) { + Emod[i][j] = Emod[j][i] = normal_coeffs_one[0]; + poiss[i][j] = poiss[j][i] = normal_coeffs_one[2]; + normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = + FOURTHIRDS*mix_stiffnessE(Emod[i][j],Emod[i][j], + poiss[i][j],poiss[i][j]); + } else { + normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = normal_coeffs_one[0]; + } + if ((normal_model_one == JKR) || (normal_model_one == DMT)) + normal_coeffs[i][j][3] = normal_coeffs[j][i][3] = normal_coeffs_one[3]; + + damping_model[i][j] = damping_model[j][i] = damping_model_one; + + tangential_model[i][j] = tangential_model[j][i] = tangential_model_one; + if (tangential_coeffs_one[0] == -1) { + tangential_coeffs[i][j][0] = tangential_coeffs[j][i][0] = + 8*mix_stiffnessG(Emod[i][j],Emod[i][j],poiss[i][j],poiss[i][j]); + } else { + tangential_coeffs[i][j][0] = tangential_coeffs[j][i][0] = + tangential_coeffs_one[0]; + } + for (int k = 1; k < 3; k++) + tangential_coeffs[i][j][k] = tangential_coeffs[j][i][k] = + tangential_coeffs_one[k]; + + roll_model[i][j] = roll_model[j][i] = roll_model_one; + if (roll_model_one != ROLL_NONE) + for (int k = 0; k < 3; k++) + roll_coeffs[i][j][k] = roll_coeffs[j][i][k] = roll_coeffs_one[k]; + + twist_model[i][j] = twist_model[j][i] = twist_model_one; + if (twist_model_one != TWIST_NONE && twist_model_one != TWIST_MARSHALL) + for (int k = 0; k < 3; k++) + twist_coeffs[i][j][k] = twist_coeffs[j][i][k] = twist_coeffs_one[k]; + + cutoff_type[i][j] = cutoff_type[j][i] = cutoff_one; + + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairGranular::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 = normal_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 = i; j <= atom->ntypes; j++) + if (normal_model[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; + } + } + for (int i = 1; i <= atom->ntypes; i++) + for (int j = i; j <= atom->ntypes; j++) + if (tangential_model[i][j] == TANGENTIAL_MINDLIN_RESCALE || + tangential_model[i][j] == TANGENTIAL_MINDLIN_RESCALE_FORCE) { + size_history += 1; + roll_history_index += 1; + twist_history_index += 1; + nondefault_history_transfer = 1; + history_transfer_factors = new int[size_history]; + for (int ii = 0; ii < size_history; ++ii) + history_transfer_factors[ii] = -1; + history_transfer_factors[3] = 1; + break; + } + + 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 and first init, create Fix to store history + // it replaces FixDummy, created in the constructor + // this is so its order in the fix list is preserved + + if (use_history && fix_history == nullptr) { + modify->replace_fix("NEIGH_HISTORY_GRANULAR_DUMMY","NEIGH_HISTORY_GRANULAR" + " all NEIGH_HISTORY " + std::to_string(size_history),1); + int ifix = modify->find_fix("NEIGH_HISTORY_GRANULAR"); + fix_history = (FixNeighHistory *) modify->fix[ifix]; + 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 = nullptr; + 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)); + onerad_dynamic[i] = radmax; + } + if (idep >= 0) { + itype = i; + double radmax = *((double *) modify->fix[idep]->extract("radius",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 (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_GRANULAR"); + 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 PairGranular::init_one(int i, int j) +{ + double cutoff=0.0; + + if (setflag[i][j] == 0) { + if ((normal_model[i][i] != normal_model[j][j]) || + (damping_model[i][i] != damping_model[j][j]) || + (tangential_model[i][i] != tangential_model[j][j]) || + (roll_model[i][i] != roll_model[j][j]) || + (twist_model[i][i] != twist_model[j][j])) { + + char str[512]; + sprintf(str,"Granular pair style functional forms are different, " + "cannot mix coefficients for types %d and %d. \n" + "This combination must be set explicitly " + "via pair_coeff command",i,j); + error->one(FLERR,str); + } + + if (normal_model[i][j] == HERTZ || normal_model[i][j] == HOOKE) + normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = + mix_geom(normal_coeffs[i][i][0], normal_coeffs[j][j][0]); + else + normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = + mix_stiffnessE(Emod[i][i], Emod[j][j], poiss[i][i], poiss[j][j]); + + normal_coeffs[i][j][1] = normal_coeffs[j][i][1] = + mix_geom(normal_coeffs[i][i][1], normal_coeffs[j][j][1]); + if ((normal_model[i][j] == JKR) || (normal_model[i][j] == 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] = tangential_coeffs[j][i][k] = + mix_geom(tangential_coeffs[i][i][k], tangential_coeffs[j][j][k]); + + if (roll_model[i][j] != 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_model[i][j] != TWIST_NONE && twist_model[i][j] != 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. + + double pulloff; + + if (cutoff_type[i][j] < 0 && cutoff_global < 0) { + if (((maxrad_dynamic[i] > 0.0) && (maxrad_dynamic[j] > 0.0)) || + ((maxrad_dynamic[i] > 0.0) && (maxrad_frozen[j] > 0.0)) || + // radius info about both i and j exist + ((maxrad_frozen[i] > 0.0) && (maxrad_dynamic[j] > 0.0))) { + cutoff = maxrad_dynamic[i]+maxrad_dynamic[j]; + pulloff = 0.0; + if (normal_model[i][j] == JKR) { + pulloff = pulloff_distance(maxrad_dynamic[i], maxrad_dynamic[j], i, j); + cutoff += pulloff; + } + + if (normal_model[i][j] == JKR) + pulloff = pulloff_distance(maxrad_frozen[i], maxrad_dynamic[j], i, j); + cutoff = MAX(cutoff, maxrad_frozen[i]+maxrad_dynamic[j]+pulloff); + + if (normal_model[i][j] == JKR) + pulloff = pulloff_distance(maxrad_dynamic[i], maxrad_frozen[j], i, j); + cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]+pulloff); + + } 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 if (cutoff_type[i][j] > 0) { + cutoff = cutoff_type[i][j]; + } else if (cutoff_global > 0) { + cutoff = cutoff_global; + } + + return cutoff; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairGranular::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_model[i][j],sizeof(int),1,fp); + fwrite(&damping_model[i][j],sizeof(int),1,fp); + fwrite(&tangential_model[i][j],sizeof(int),1,fp); + fwrite(&roll_model[i][j],sizeof(int),1,fp); + fwrite(&twist_model[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(&cutoff_type[i][j],sizeof(double),1,fp); + } + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairGranular::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) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error); + MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (setflag[i][j]) { + if (me == 0) { + utils::sfread(FLERR,&normal_model[i][j],sizeof(int),1,fp,nullptr,error); + utils::sfread(FLERR,&damping_model[i][j],sizeof(int),1,fp,nullptr,error); + utils::sfread(FLERR,&tangential_model[i][j],sizeof(int),1,fp,nullptr,error); + utils::sfread(FLERR,&roll_model[i][j],sizeof(int),1,fp,nullptr,error); + utils::sfread(FLERR,&twist_model[i][j],sizeof(int),1,fp,nullptr,error); + utils::sfread(FLERR,normal_coeffs[i][j],sizeof(double),4,fp,nullptr,error); + utils::sfread(FLERR,tangential_coeffs[i][j],sizeof(double),3,fp,nullptr,error); + utils::sfread(FLERR,roll_coeffs[i][j],sizeof(double),3,fp,nullptr,error); + utils::sfread(FLERR,twist_coeffs[i][j],sizeof(double),3,fp,nullptr,error); + utils::sfread(FLERR,&cutoff_type[i][j],sizeof(double),1,fp,nullptr,error); + } + MPI_Bcast(&normal_model[i][j],1,MPI_INT,0,world); + MPI_Bcast(&damping_model[i][j],1,MPI_INT,0,world); + MPI_Bcast(&tangential_model[i][j],1,MPI_INT,0,world); + MPI_Bcast(&roll_model[i][j],1,MPI_INT,0,world); + MPI_Bcast(&twist_model[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(&cutoff_type[i][j],1,MPI_DOUBLE,0,world); + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +void PairGranular::reset_dt() +{ + dt = update->dt; +} + +/* ---------------------------------------------------------------------- */ + +double PairGranular::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,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; + double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3; + + double knfac, damp_normal, damp_normal_prefactor; + double k_tangential, damp_tangential; + double Fne, Ft, Fdamp, Fntot, Fncrit, 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; + + // rolling + double k_roll, damp_roll; + double rollmag; + double fr, fr1, fr2, fr3; + + // twisting + double k_twist, damp_twist, mu_twist; + double signtwist, magtwist, magtortwist, Mtcrit; + + double shrmag; + int jnum; + int *jlist; + double *history,*allhistory; + + int nall = atom->nlocal + atom->nghost; + if ((i >= nall) || (j >= nall)) + error->all(FLERR,"Not enough atoms for pair granular single function"); + + double *radius = atom->radius; + radi = radius[i]; + radj = radius[j]; + radsum = radi + radj; + Reff = (radsum > 0.0) ? radi*radj/radsum : 0.0; + + bool touchflag; + E = normal_coeffs[itype][jtype][0]; + if (normal_model[itype][jtype] == JKR) { + E *= THREEQUARTERS; + R2 = Reff*Reff; + coh = normal_coeffs[itype][jtype][3]; + a = cbrt(9.0*MY_PI*coh*R2/(4*E)); + delta_pulloff = a*a/Reff - 2*sqrt(MY_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; + + // 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 *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; + + delta = radsum - r; + dR = delta*Reff; + if (normal_model[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; + // in case sqrt(0) < 0 due to precision issues + sqrt1 = MAX(0, t0*(t1+2*t2)); + t4 = cbrt(t1+t2+THREEROOT3*MY_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*MY_PI*R2/(E*t6)); + a = INVROOT6*(t6 + sqrt(sqrt3)); + a2 = a*a; + knfac = normal_coeffs[itype][jtype][0]*a; + Fne = knfac*a2/Reff - MY_2PI*a2*sqrt(4*coh*E/(MY_PI*a)); + } else { + knfac = E; + Fne = knfac*delta; + a = sqrt(dR); + if (normal_model[itype][jtype] != HOOKE) { + Fne *= a; + knfac *= a; + } + if (normal_model[itype][jtype] == DMT) + Fne -= 4*MY_PI*normal_coeffs[itype][jtype][3]*Reff; + } + + if (damping_model[itype][jtype] == VELOCITY) { + damp_normal = 1; + } else if (damping_model[itype][jtype] == MASS_VELOCITY) { + damp_normal = meff; + } else if (damping_model[itype][jtype] == VISCOELASTIC) { + damp_normal = a*meff; + } else if (damping_model[itype][jtype] == TSUJI) { + damp_normal = sqrt(meff*knfac); + } else damp_normal = 0.0; + + damp_normal_prefactor = normal_coeffs[itype][jtype][1]*damp_normal; + Fdamp = -damp_normal_prefactor*vnnr; + + Fntot = Fne + Fdamp; + + jnum = list->numneigh[i]; + jlist = list->firstneigh[i]; + + if (use_history) { + if ((fix_history == nullptr) || (fix_history->firstvalue == nullptr)) + error->one(FLERR,"Pair granular single computation needs 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 + //**************************************** + + // For linear, mindlin, mindlin_rescale: + // history = cumulative tangential displacement + // + // For mindlin/force, mindlin_rescale/force: + // history = cumulative tangential elastic force + + // 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_model[itype][jtype] == JKR) { + F_pulloff = 3*MY_PI*coh*Reff; + Fncrit = fabs(Fne + 2*F_pulloff); + } else if (normal_model[itype][jtype] == DMT) { + F_pulloff = 4*MY_PI*coh*Reff; + Fncrit = fabs(Fne + 2*F_pulloff); + } else { + Fncrit = fabs(Fntot); + } + Fscrit = tangential_coeffs[itype][jtype][2] * Fncrit; + + //------------------------------ + // tangential forces + //------------------------------ + k_tangential = tangential_coeffs[itype][jtype][0]; + damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal_prefactor; + + if (tangential_history) { + if (tangential_model[itype][jtype] != TANGENTIAL_HISTORY) { + k_tangential *= a; + } + + shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + + history[2]*history[2]); + + // tangential forces = history + tangential velocity damping + if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY || + tangential_model[itype][jtype] == TANGENTIAL_MINDLIN || + tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) { + fs1 = -k_tangential*history[0] - damp_tangential*vtr1; + fs2 = -k_tangential*history[1] - damp_tangential*vtr2; + fs3 = -k_tangential*history[2] - damp_tangential*vtr3; + } else { + fs1 = history[0] - damp_tangential*vtr1; + fs2 = history[1] - damp_tangential*vtr2; + fs3 = history[2] - damp_tangential*vtr3; + } + + // rescale frictional forces if needed + fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); + 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; + } + + // classic pair gran/hooke (no history) + } else { + fs = damp_tangential*vrel; + if (vrel != 0.0) Ft = MIN(Fscrit,fs) / vrel; + else Ft = 0.0; + fs1 = -Ft*vtr1; + fs2 = -Ft*vtr2; + fs3 = -Ft*vtr3; + fs = Ft*vrel; + } + + //**************************************** + // rolling resistance + //**************************************** + + if ((roll_model[itype][jtype] != ROLL_NONE) + || (twist_model[itype][jtype] != TWIST_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; + + 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]); + + 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] * Fncrit; + + fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3); + if (fr > Frcrit) { + if (rollmag != 0.0) { + fr1 *= Frcrit/fr; + fr2 *= Frcrit/fr; + fr3 *= Frcrit/fr; + fr *= Frcrit/fr; + } else fr1 = fr2 = fr3 = fr = 0.0; + } + } else fr1 = fr2 = fr3 = fr = 0.0; + + //**************************************** + // twisting torque, including history effects + //**************************************** + + if (twist_model[itype][jtype] != TWIST_NONE) { + // omega_T (eq 29 of Marshall) + magtwist = relrot1*nx + relrot2*ny + relrot3*nz; + if (twist_model[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*tangential_coeffs[itype][jtype][2];; + } else { + k_twist = twist_coeffs[itype][jtype][0]; + damp_twist = twist_coeffs[itype][jtype][1]; + mu_twist = twist_coeffs[itype][jtype][2]; + } + // M_t torque (eq 30) + magtortwist = -k_twist*history[twist_history_index] - damp_twist*magtwist; + signtwist = (magtwist > 0) - (magtwist < 0); + Mtcrit = mu_twist*Fncrit; // critical torque (eq 44) + if (fabs(magtortwist) > Mtcrit) { + magtortwist = -Mtcrit * signtwist; // eq 34 + } else magtortwist = 0.0; + } else magtortwist = 0.0; + + // set force and return no energy + + fforce = Fntot*rinv; + + // 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; + svector[9] = delx; + svector[10] = dely; + svector[11] = delz; + return 0.0; +} + +/* ---------------------------------------------------------------------- */ + +int PairGranular::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 PairGranular::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 PairGranular::memory_usage() +{ + double bytes = (double)nmax * sizeof(double); + return bytes; +} + +/* ---------------------------------------------------------------------- + mixing of Young's modulus (E) +------------------------------------------------------------------------- */ + +double PairGranular::mix_stiffnessE(double Eii, double Ejj, + double poisii, double poisjj) +{ + return 1/((1-poisii*poisii)/Eii+(1-poisjj*poisjj)/Ejj); +} + +/* ---------------------------------------------------------------------- + mixing of shear modulus (G) +------------------------------------------------------------------------ */ + +double PairGranular::mix_stiffnessG(double Eii, double Ejj, + double poisii, double poisjj) +{ + return 1/((2*(2-poisii)*(1+poisii)/Eii) + (2*(2-poisjj)*(1+poisjj)/Ejj)); +} + +/* ---------------------------------------------------------------------- + mixing of everything else +------------------------------------------------------------------------- */ + +double PairGranular::mix_geom(double valii, double valjj) +{ + return sqrt(valii*valjj); +} + + +/* ---------------------------------------------------------------------- + compute pull-off distance (beyond contact) for a given radius and atom type +------------------------------------------------------------------------- */ + +double PairGranular::pulloff_distance(double radi, double radj, + int itype, int jtype) +{ + double E, coh, a, Reff; + Reff = radi*radj/(radi+radj); + if (Reff <= 0) return 0; + coh = normal_coeffs[itype][jtype][3]; + E = normal_coeffs[itype][jtype][0]*THREEQUARTERS; + a = cbrt(9*MY_PI*coh*Reff*Reff/(4*E)); + return a*a/Reff - 2*sqrt(MY_PI*coh*a/E); +} + +/* ---------------------------------------------------------------------- + transfer history during fix/neigh/history exchange + only needed if any history entries i-j are not just negative of j-i entries +------------------------------------------------------------------------- */ + +void PairGranular::transfer_history(double* source, double* target) +{ + for (int i = 0; i < size_history; i++) + target[i] = history_transfer_factors[i]*source[i]; +} From 773ec40d3a83008e0e9bcfc0d51443ba9b6e2dd0 Mon Sep 17 00:00:00 2001 From: Joel Clemmer Date: Mon, 5 Apr 2021 11:37:35 -0600 Subject: [PATCH 07/15] Misc small fixes --- src/GRANULAR/fix_wall_gran.cpp | 11 +- src/GRANULAR/pair_gran_hertz_history.cpp | 4 +- src/GRANULAR/pair_gran_hooke.cpp | 4 +- src/GRANULAR/pair_granular.cpp | 2 +- src/GRANULAR/pair_granular_corrected.cpp | 1821 ---------------------- 5 files changed, 10 insertions(+), 1832 deletions(-) delete mode 100644 src/GRANULAR/pair_granular_corrected.cpp diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index 1d27e21162..7ce8e49008 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -376,7 +376,6 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : wiggle = 0; wshear = 0; peratom_flag = 0; - limit_damping = 0; while (iarg < narg) { if (strcmp(arg[iarg],"wiggle") == 0) { @@ -792,7 +791,7 @@ void FixWallGran::hooke(double rsq, double dx, double dy, double dz, damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radius-r)*rinv - damp; - if(limit_damping and ccel < 0.0) ccel = 0.0; + if (limit_damping and ccel < 0.0) ccel = 0.0; // relative velocities @@ -885,7 +884,7 @@ void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz, damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radius-r)*rinv - damp; - if(limit_damping and ccel < 0.0) ccel = 0.0; + if (limit_damping and ccel < 0.0) ccel = 0.0; // relative velocities @@ -1017,7 +1016,7 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz, if (rwall == 0.0) polyhertz = sqrt((radius-r)*radius); else polyhertz = sqrt((radius-r)*radius*rwall/(rwall+radius)); ccel *= polyhertz; - if(limit_damping and ccel < 0.0) ccel = 0.0; + if (limit_damping and ccel < 0.0) ccel = 0.0; // relative velocities @@ -1303,8 +1302,8 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz, history[thist2] -= rsht*nz; // also rescale to preserve magnitude - prjmag = sqrt(history[0]*history[0] + history[1]*history[1] + - history[2]*history[2]); + prjmag = sqrt(history[thist0]*history[thist0] + + history[thist1]*history[thist1] + history[thist2]*history[thist2]); if (prjmag > 0) scalefac = shrmag/prjmag; else scalefac = 0; history[thist0] *= scalefac; diff --git a/src/GRANULAR/pair_gran_hertz_history.cpp b/src/GRANULAR/pair_gran_hertz_history.cpp index 40c6c3a41d..e21ec727d6 100644 --- a/src/GRANULAR/pair_gran_hertz_history.cpp +++ b/src/GRANULAR/pair_gran_hertz_history.cpp @@ -181,7 +181,7 @@ void PairGranHertzHistory::compute(int eflag, int vflag) ccel = kn*(radsum-r)*rinv - damp; polyhertz = sqrt((radsum-r)*radi*radj / radsum); ccel *= polyhertz; - if(limit_damping and ccel < 0.0) ccel = 0.0; + if (limit_damping and ccel < 0.0) ccel = 0.0; // relative velocities @@ -389,7 +389,7 @@ double PairGranHertzHistory::single(int i, int j, int /*itype*/, int /*jtype*/, ccel = kn*(radsum-r)*rinv - damp; polyhertz = sqrt((radsum-r)*radi*radj / radsum); ccel *= polyhertz; - if(limit_damping and ccel < 0.0) ccel = 0.0; + if (limit_damping and ccel < 0.0) ccel = 0.0; // relative velocities diff --git a/src/GRANULAR/pair_gran_hooke.cpp b/src/GRANULAR/pair_gran_hooke.cpp index 0e3ffc7293..c28bbd007f 100644 --- a/src/GRANULAR/pair_gran_hooke.cpp +++ b/src/GRANULAR/pair_gran_hooke.cpp @@ -158,7 +158,7 @@ void PairGranHooke::compute(int eflag, int vflag) damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; - if(limit_damping and ccel < 0.0) ccel = 0.0; + if (limit_damping and ccel < 0.0) ccel = 0.0; // relative velocities @@ -300,7 +300,7 @@ double PairGranHooke::single(int i, int j, int /*itype*/, int /*jtype*/, double damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; - if(limit_damping and ccel < 0.0) ccel = 0.0; + if (limit_damping and ccel < 0.0) ccel = 0.0; // relative velocities diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 8bf9b0c4ed..ac98a6e27c 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -1547,7 +1547,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, Fdamp = -damp_normal_prefactor*vnnr; Fntot = Fne + Fdamp; - if(limit_damping[itype][jtype] and Fntot < 0.0) Fntot = 0.0; + if (limit_damping[itype][jtype] and Fntot < 0.0) Fntot = 0.0; jnum = list->numneigh[i]; jlist = list->firstneigh[i]; diff --git a/src/GRANULAR/pair_granular_corrected.cpp b/src/GRANULAR/pair_granular_corrected.cpp deleted file mode 100644 index bcf262aa2c..0000000000 --- a/src/GRANULAR/pair_granular_corrected.cpp +++ /dev/null @@ -1,1821 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://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 "pair_granular.h" - -#include -#include - -#include "atom.h" -#include "force.h" -#include "update.h" -#include "modify.h" -#include "fix.h" -#include "fix_dummy.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 "math_special.h" - - -using namespace LAMMPS_NS; -using namespace MathConst; -using namespace MathSpecial; - -#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 4.0/3.0 // 4/3 -#define THREEQUARTERS 0.75 // 3/4 - -#define EPSILON 1e-10 - -enum {HOOKE, HERTZ, HERTZ_MATERIAL, DMT, JKR}; -enum {VELOCITY, MASS_VELOCITY, VISCOELASTIC, TSUJI}; -enum {TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY, - TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE, - TANGENTIAL_MINDLIN_FORCE, TANGENTIAL_MINDLIN_RESCALE_FORCE}; -enum {TWIST_NONE, TWIST_SDS, TWIST_MARSHALL}; -enum {ROLL_NONE, ROLL_SDS}; - -/* ---------------------------------------------------------------------- */ - -PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp) -{ - single_enable = 1; - no_virial_fdotr_compute = 1; - centroidstressflag = CENTROID_NOTAVAIL; - - single_extra = 12; - svector = new double[single_extra]; - - neighprev = 0; - - nmax = 0; - mass_rigid = nullptr; - - onerad_dynamic = nullptr; - onerad_frozen = nullptr; - maxrad_dynamic = nullptr; - maxrad_frozen = nullptr; - - history_transfer_factors = nullptr; - - 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; - - // create dummy fix as placeholder for FixNeighHistory - // this is so final order of Modify:fix will conform to input script - - fix_history = nullptr; - modify->add_fix("NEIGH_HISTORY_GRANULAR_DUMMY all DUMMY"); - fix_dummy = (FixDummy *) modify->fix[modify->nfix-1]; -} - -/* ---------------------------------------------------------------------- */ - -PairGranular::~PairGranular() -{ - delete [] svector; - - if (!fix_history) modify->delete_fix("NEIGH_HISTORY_GRANULAR_DUMMY"); - else modify->delete_fix("NEIGH_HISTORY_GRANULAR"); - - if (allocated) { - memory->destroy(setflag); - memory->destroy(cutsq); - memory->destroy(cutoff_type); - - memory->destroy(normal_coeffs); - memory->destroy(tangential_coeffs); - memory->destroy(roll_coeffs); - memory->destroy(twist_coeffs); - - memory->destroy(Emod); - memory->destroy(poiss); - - memory->destroy(normal_model); - memory->destroy(damping_model); - memory->destroy(tangential_model); - memory->destroy(roll_model); - memory->destroy(twist_model); - - delete [] onerad_dynamic; - delete [] onerad_frozen; - delete [] maxrad_dynamic; - delete [] maxrad_frozen; - } - - memory->destroy(mass_rigid); -} - -/* ---------------------------------------------------------------------- */ - -void PairGranular::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; - double Reff, delta, dR, dR2, dist_to_contact; - - double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; - double wr1,wr2,wr3; - double vtr1,vtr2,vtr3,vrel; - - double knfac, damp_normal=0.0, damp_normal_prefactor; - double k_tangential, damp_tangential; - double Fne, Ft, Fdamp, Fntot, Fncrit, Fscrit, Frcrit; - double fs, fs1, fs2, fs3, tor1, tor2, tor3; - - double mi,mj,meff; - double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3; - - // 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; - - // rolling - double k_roll, damp_roll; - double 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,prjmag; - bool frameupdate; - int *ilist,*jlist,*numneigh,**firstneigh; - int *touch,**firsttouch; - double *history,*allhistory,**firsthistory; - - bool touchflag = false; - const bool historyupdate = (update->setupflag) ? false : true; - - ev_init(eflag,vflag); - - // 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; - - inum = list->inum; - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->firstneigh; - if (use_history) { - 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]; - if (use_history) { - 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/radsum; - touchflag = false; - - if (normal_model[itype][jtype] == JKR) { - E *= THREEQUARTERS; - if (touch[jj]) { - R2 = Reff*Reff; - coh = normal_coeffs[itype][jtype][3]; - a = cbrt(9.0*MY_PI*coh*R2/(4*E)); - delta_pulloff = a*a/Reff - 2*sqrt(MY_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 - if (use_history) { - 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_model[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; - // in case sqrt(0) < 0 due to precision issues - sqrt1 = MAX(0, t0*(t1+2*t2)); - t4 = cbrt(t1+t2+THREEROOT3*MY_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*MY_PI*R2/(E*t6)); - a = INVROOT6*(t6 + sqrt(sqrt3)); - a2 = a*a; - knfac = normal_coeffs[itype][jtype][0]*a; - Fne = knfac*a2/Reff - MY_2PI*a2*sqrt(4*coh*E/(MY_PI*a)); - } else { - knfac = E; // Hooke - Fne = knfac*delta; - a = sqrt(dR); - if (normal_model[itype][jtype] != HOOKE) { - Fne *= a; - knfac *= a; - } - if (normal_model[itype][jtype] == DMT) - Fne -= 4*MY_PI*normal_coeffs[itype][jtype][3]*Reff; - } - - // NOTE: consider restricting Hooke to only have - // 'velocity' as an option for damping? - - if (damping_model[itype][jtype] == VELOCITY) { - damp_normal = 1; - } else if (damping_model[itype][jtype] == MASS_VELOCITY) { - damp_normal = meff; - } else if (damping_model[itype][jtype] == VISCOELASTIC) { - damp_normal = a*meff; - } else if (damping_model[itype][jtype] == TSUJI) { - damp_normal = sqrt(meff*knfac); - } - - damp_normal_prefactor = normal_coeffs[itype][jtype][1]*damp_normal; - Fdamp = -damp_normal_prefactor*vnnr; - - Fntot = Fne + Fdamp; - - //**************************************** - // tangential force, including history effects - //**************************************** - - // For linear, mindlin, mindlin_rescale: - // history = cumulative tangential displacement - // - // For mindlin/force, mindlin_rescale/force: - // history = cumulative tangential elastic force - - // 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_model[itype][jtype] == JKR) { - F_pulloff = 3*MY_PI*coh*Reff; - Fncrit = fabs(Fne + 2*F_pulloff); - } else if (normal_model[itype][jtype] == DMT) { - F_pulloff = 4*MY_PI*coh*Reff; - Fncrit = fabs(Fne + 2*F_pulloff); - } else { - Fncrit = fabs(Fntot); - } - Fscrit = tangential_coeffs[itype][jtype][2] * Fncrit; - - //------------------------------ - // tangential forces - //------------------------------ - k_tangential = tangential_coeffs[itype][jtype][0]; - damp_tangential = tangential_coeffs[itype][jtype][1] * - damp_normal_prefactor; - - if (tangential_history) { - if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN || - tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_FORCE) { - k_tangential *= a; - } else if (tangential_model[itype][jtype] == - TANGENTIAL_MINDLIN_RESCALE || - tangential_model[itype][jtype] == - TANGENTIAL_MINDLIN_RESCALE_FORCE) { - k_tangential *= a; - // on unloading, rescale the shear displacements/force - if (a < history[3]) { - double factor = a/history[3]; - history[0] *= factor; - history[1] *= factor; - history[2] *= factor; - } - } - // rotate and update displacements / force. - // 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 (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_FORCE || - tangential_model[itype][jtype] == - TANGENTIAL_MINDLIN_RESCALE_FORCE) - frameupdate = fabs(rsht) > EPSILON*Fscrit; - else - frameupdate = fabs(rsht)*k_tangential > EPSILON*Fscrit; - if (frameupdate) { - shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + - history[2]*history[2]); - // projection - history[0] -= rsht*nx; - history[1] -= rsht*ny; - history[2] -= rsht*nz; - - // also rescale to preserve magnitude - prjmag = sqrt(history[0]*history[0] + history[1]*history[1] + - history[2]*history[2]); - if (prjmag > 0) scalefac = shrmag/prjmag; - else scalefac = 0; - history[0] *= scalefac; - history[1] *= scalefac; - history[2] *= scalefac; - } - // update history - if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY || - tangential_model[itype][jtype] == TANGENTIAL_MINDLIN || - tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) { - // tangential displacement - history[0] += vtr1*dt; - history[1] += vtr2*dt; - history[2] += vtr3*dt; - } else { - // tangential force - // see e.g. eq. 18 of Thornton et al, Pow. Tech. 2013, v223,p30-46 - history[0] -= k_tangential*vtr1*dt; - history[1] -= k_tangential*vtr2*dt; - history[2] -= k_tangential*vtr3*dt; - } - if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE || - tangential_model[itype][jtype] == - TANGENTIAL_MINDLIN_RESCALE_FORCE) - history[3] = a; - } - - // tangential forces = history + tangential velocity damping - if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY || - tangential_model[itype][jtype] == TANGENTIAL_MINDLIN || - tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) { - fs1 = -k_tangential*history[0] - damp_tangential*vtr1; - fs2 = -k_tangential*history[1] - damp_tangential*vtr2; - fs3 = -k_tangential*history[2] - damp_tangential*vtr3; - } else { - fs1 = history[0] - damp_tangential*vtr1; - fs2 = history[1] - damp_tangential*vtr2; - fs3 = history[2] - damp_tangential*vtr3; - } - - // rescale frictional displacements and forces if needed - fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); - if (fs > Fscrit) { - shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + - history[2]*history[2]); - if (shrmag != 0.0) { - if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY || - tangential_model[itype][jtype] == TANGENTIAL_MINDLIN || - tangential_model[itype][jtype] == - TANGENTIAL_MINDLIN_RESCALE) { - 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); - } else { - history[0] = Fscrit*fs1/fs + damp_tangential*vtr1; - history[1] = Fscrit*fs2/fs + damp_tangential*vtr2; - history[2] = 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 = damp_tangential*vrel; - if (vrel != 0.0) Ft = MIN(Fscrit,fs) / vrel; - else Ft = 0.0; - fs1 = -Ft*vtr1; - fs2 = -Ft*vtr2; - fs3 = -Ft*vtr3; - } - - if (roll_model[itype][jtype] != ROLL_NONE || - twist_model[itype][jtype] != TWIST_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) - // - 0.5*((radj-radi)/radsum)*vtr1; - // - 0.5*((radj-radi)/radsum)*vtr2; - // - 0.5*((radj-radi)/radsum)*vtr3; - } - //**************************************** - // rolling resistance - //**************************************** - - if (roll_model[itype][jtype] != ROLL_NONE) { - vrl1 = Reff*(relrot2*nz - relrot3*ny); - vrl2 = Reff*(relrot3*nx - relrot1*nz); - vrl3 = Reff*(relrot1*ny - relrot2*nx); - - int rhist0 = roll_history_index; - int rhist1 = rhist0 + 1; - int rhist2 = rhist1 + 1; - - k_roll = roll_coeffs[itype][jtype][0]; - damp_roll = roll_coeffs[itype][jtype][1]; - Frcrit = roll_coeffs[itype][jtype][2] * Fncrit; - - if (historyupdate) { - rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz; - frameupdate = fabs(rolldotn)*k_roll > EPSILON*Frcrit; - if (frameupdate) { // rotate into tangential plane - rollmag = sqrt(history[rhist0]*history[rhist0] + - history[rhist1]*history[rhist1] + - history[rhist2]*history[rhist2]); - // projection - history[rhist0] -= rolldotn*nx; - history[rhist1] -= rolldotn*ny; - history[rhist2] -= rolldotn*nz; - // also rescale to preserve magnitude - prjmag = sqrt(history[rhist0]*history[rhist0] + - history[rhist1]*history[rhist1] + - history[rhist2]*history[rhist2]); - if (prjmag > 0) scalefac = rollmag/prjmag; - else scalefac = 0; - history[rhist0] *= scalefac; - history[rhist1] *= scalefac; - history[rhist2] *= scalefac; - } - history[rhist0] += vrl1*dt; - history[rhist1] += vrl2*dt; - history[rhist2] += vrl3*dt; - } - - 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 - - fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3); - if (fr > Frcrit) { - rollmag = sqrt(history[rhist0]*history[rhist0] + - history[rhist1]*history[rhist1] + - history[rhist2]*history[rhist2]); - 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; - } - } - - //**************************************** - // twisting torque, including history effects - //**************************************** - - if (twist_model[itype][jtype] != TWIST_NONE) { - // omega_T (eq 29 of Marshall) - magtwist = relrot1*nx + relrot2*ny + relrot3*nz; - if (twist_model[itype][jtype] == TWIST_MARSHALL) { - k_twist = 0.5*k_tangential*a*a;; // eq 32 of Marshall paper - damp_twist = 0.5*damp_tangential*a*a; - mu_twist = TWOTHIRDS*a*tangential_coeffs[itype][jtype][2]; - } else { - k_twist = twist_coeffs[itype][jtype][0]; - damp_twist = twist_coeffs[itype][jtype][1]; - mu_twist = twist_coeffs[itype][jtype][2]; - } - 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 = mu_twist*Fncrit; // 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 - } - } - - // 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; - - dist_to_contact = radi-0.5*delta; - torque[i][0] -= dist_to_contact*tor1; - torque[i][1] -= dist_to_contact*tor2; - torque[i][2] -= dist_to_contact*tor3; - - if (twist_model[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_model[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; - - dist_to_contact = radj-0.5*delta; - torque[j][0] -= dist_to_contact*tor1; - torque[j][1] -= dist_to_contact*tor2; - torque[j][2] -= dist_to_contact*tor3; - - if (twist_model[itype][jtype] != TWIST_NONE) { - torque[j][0] -= tortwist1; - torque[j][1] -= tortwist2; - torque[j][2] -= tortwist3; - } - if (roll_model[itype][jtype] != ROLL_NONE) { - torque[j][0] -= torroll1; - torque[j][1] -= torroll2; - torque[j][2] -= torroll3; - } - } - if (evflag) ev_tally_xyz(i,j,nlocal,force->newton_pair, - 0.0,0.0,fx,fy,fz,delx,dely,delz); - } - } - } -} - -/* ---------------------------------------------------------------------- - allocate all arrays -------------------------------------------------------------------------- */ - -void PairGranular::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(cutoff_type,n+1,n+1,"pair:cutoff_type"); - 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(Emod,n+1,n+1,"pair:Emod"); - memory->create(poiss,n+1,n+1,"pair:poiss"); - - memory->create(normal_model,n+1,n+1,"pair:normal_model"); - memory->create(damping_model,n+1,n+1,"pair:damping_model"); - memory->create(tangential_model,n+1,n+1,"pair:tangential_model"); - memory->create(roll_model,n+1,n+1,"pair:roll_model"); - memory->create(twist_model,n+1,n+1,"pair:twist_model"); - - 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 PairGranular::settings(int narg, char **arg) -{ - if (narg == 1) { - cutoff_global = utils::numeric(FLERR,arg[0],false,lmp); - } else { - cutoff_global = -1; // will be set based on particle sizes, model choice - } - - normal_history = tangential_history = 0; - roll_history = twist_history = 0; -} - -/* ---------------------------------------------------------------------- - set coeffs for one or more type pairs -------------------------------------------------------------------------- */ - -void PairGranular::coeff(int narg, char **arg) -{ - int normal_model_one, damping_model_one; - int tangential_model_one, roll_model_one, twist_model_one; - - double normal_coeffs_one[4]; - double tangential_coeffs_one[4]; - double roll_coeffs_one[4]; - double twist_coeffs_one[4]; - - double cutoff_one = -1; - - if (narg < 2) - error->all(FLERR,"Incorrect args for pair coefficients"); - - if (!allocated) allocate(); - - int ilo,ihi,jlo,jhi; - utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); - utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); - - //Defaults - normal_model_one = tangential_model_one = -1; - roll_model_one = ROLL_NONE; - twist_model_one = TWIST_NONE; - damping_model_one = 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_model_one = HOOKE; - normal_coeffs_one[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp); // kn - normal_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp); // damping - iarg += 3; - } else if (strcmp(arg[iarg], "hertz") == 0) { - if (iarg + 2 >= narg) - error->all(FLERR,"Illegal pair_coeff command, " - "not enough parameters provided for Hertz option"); - normal_model_one = HERTZ; - normal_coeffs_one[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp); // kn - normal_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp); // damping - iarg += 3; - } else if (strcmp(arg[iarg], "hertz/material") == 0) { - if (iarg + 3 >= narg) - error->all(FLERR,"Illegal pair_coeff command, " - "not enough parameters provided for Hertz/material option"); - normal_model_one = HERTZ_MATERIAL; - normal_coeffs_one[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp); // E - normal_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp); // damping - normal_coeffs_one[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp); // Poisson's ratio - iarg += 4; - } 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_model_one = DMT; - normal_coeffs_one[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp); // E - normal_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp); // damping - normal_coeffs_one[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp); // Poisson's ratio - normal_coeffs_one[3] = utils::numeric(FLERR,arg[iarg+4],false,lmp); // 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_model_one = JKR; - normal_coeffs_one[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp); // E - normal_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp); // damping - normal_coeffs_one[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp); // Poisson's ratio - normal_coeffs_one[3] = utils::numeric(FLERR,arg[iarg+4],false,lmp); // cohesion - iarg += 5; - } else if (strcmp(arg[iarg], "damping") == 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_model_one = VELOCITY; - iarg += 1; - } else if (strcmp(arg[iarg+1], "mass_velocity") == 0) { - damping_model_one = MASS_VELOCITY; - iarg += 1; - } else if (strcmp(arg[iarg+1], "viscoelastic") == 0) { - damping_model_one = VISCOELASTIC; - iarg += 1; - } else if (strcmp(arg[iarg+1], "tsuji") == 0) { - damping_model_one = TSUJI; - iarg += 1; - } else error->all(FLERR, "Illegal pair_coeff command, " - "unrecognized damping model"); - iarg += 1; - } else if (strcmp(arg[iarg], "tangential") == 0) { - if (iarg + 1 >= narg) - error->all(FLERR,"Illegal pair_coeff command, must specify " - "tangential model after tangential keyword"); - if (strcmp(arg[iarg+1], "linear_nohistory") == 0) { - if (iarg + 3 >= narg) - error->all(FLERR,"Illegal pair_coeff command, " - "not enough parameters provided for tangential model"); - tangential_model_one = TANGENTIAL_NOHISTORY; - tangential_coeffs_one[0] = 0; - // gammat and friction coeff - tangential_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp); - tangential_coeffs_one[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp); - iarg += 4; - } else if ((strcmp(arg[iarg+1], "linear_history") == 0) || - (strcmp(arg[iarg+1], "mindlin") == 0) || - (strcmp(arg[iarg+1], "mindlin_rescale") == 0) || - (strcmp(arg[iarg+1], "mindlin/force") == 0) || - (strcmp(arg[iarg+1], "mindlin_rescale/force") == 0)) { - if (iarg + 4 >= narg) - error->all(FLERR,"Illegal pair_coeff command, " - "not enough parameters provided for tangential model"); - if (strcmp(arg[iarg+1], "linear_history") == 0) - tangential_model_one = TANGENTIAL_HISTORY; - else if (strcmp(arg[iarg+1], "mindlin") == 0) - tangential_model_one = TANGENTIAL_MINDLIN; - else if (strcmp(arg[iarg+1], "mindlin_rescale") == 0) - tangential_model_one = TANGENTIAL_MINDLIN_RESCALE; - else if (strcmp(arg[iarg+1], "mindlin/force") == 0) - tangential_model_one = TANGENTIAL_MINDLIN_FORCE; - else if (strcmp(arg[iarg+1], "mindlin_rescale/force") == 0) - tangential_model_one = TANGENTIAL_MINDLIN_RESCALE_FORCE; - tangential_history = 1; - if ((tangential_model_one == TANGENTIAL_MINDLIN || - tangential_model_one == TANGENTIAL_MINDLIN_RESCALE || - tangential_model_one == TANGENTIAL_MINDLIN_FORCE || - tangential_model_one == TANGENTIAL_MINDLIN_RESCALE_FORCE) && - (strcmp(arg[iarg+2], "NULL") == 0)) { - if (normal_model_one == HERTZ || normal_model_one == HOOKE) { - error->all(FLERR, "NULL setting for Mindlin tangential " - "stiffness requires a normal contact model that " - "specifies material properties"); - } - tangential_coeffs_one[0] = -1; - } else { - tangential_coeffs_one[0] = utils::numeric(FLERR,arg[iarg+2],false,lmp); // kt - } - // gammat and friction coeff - tangential_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+3],false,lmp); - tangential_coeffs_one[2] = utils::numeric(FLERR,arg[iarg+4],false,lmp); - iarg += 5; - } else { - error->all(FLERR, "Illegal pair_coeff command, " - "tangential model not recognized"); - } - } 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_model_one = ROLL_NONE; - iarg += 2; - } else if (strcmp(arg[iarg+1], "sds") == 0) { - if (iarg + 4 >= narg) - error->all(FLERR,"Illegal pair_coeff command, " - "not enough parameters provided for rolling model"); - roll_model_one = ROLL_SDS; - roll_history = 1; - // kR and gammaR and rolling friction coeff - roll_coeffs_one[0] = utils::numeric(FLERR,arg[iarg+2],false,lmp); - roll_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+3],false,lmp); - roll_coeffs_one[2] = utils::numeric(FLERR,arg[iarg+4],false,lmp); - iarg += 5; - } else { - error->all(FLERR, "Illegal pair_coeff command, " - "rolling friction model not recognized"); - } - } 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_model_one = TWIST_NONE; - iarg += 2; - } else if (strcmp(arg[iarg+1], "marshall") == 0) { - twist_model_one = TWIST_MARSHALL; - twist_history = 1; - iarg += 2; - } else if (strcmp(arg[iarg+1], "sds") == 0) { - if (iarg + 4 >= narg) - error->all(FLERR,"Illegal pair_coeff command, " - "not enough parameters provided for twist model"); - twist_model_one = TWIST_SDS; - twist_history = 1; - // kt and gammat and friction coeff - twist_coeffs_one[0] = utils::numeric(FLERR,arg[iarg+2],false,lmp); - twist_coeffs_one[1] = utils::numeric(FLERR,arg[iarg+3],false,lmp); - twist_coeffs_one[2] = utils::numeric(FLERR,arg[iarg+4],false,lmp); - iarg += 5; - } else { - error->all(FLERR, "Illegal pair_coeff command, " - "twisting friction model not recognized"); - } - } else if (strcmp(arg[iarg], "cutoff") == 0) { - if (iarg + 1 >= narg) - error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); - cutoff_one = utils::numeric(FLERR,arg[iarg+1],false,lmp); - iarg += 2; - } else error->all(FLERR, "Illegal pair coeff command"); - } - - // error not to specify normal or tangential model - if ((normal_model_one < 0) || (tangential_model_one < 0)) - error->all(FLERR, "Illegal pair coeff command, " - "must specify normal or tangential contact model"); - - int count = 0; - double damp; - if (damping_model_one == TSUJI) { - double cor; - cor = normal_coeffs_one[1]; - damp = 1.2728-4.2783*cor+11.087*square(cor)-22.348*cube(cor)+ - 27.467*powint(cor,4)-18.022*powint(cor,5)+4.8218*powint(cor,6); - } else damp = normal_coeffs_one[1]; - - for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { - normal_model[i][j] = normal_model[j][i] = normal_model_one; - normal_coeffs[i][j][1] = normal_coeffs[j][i][1] = damp; - if (normal_model_one != HERTZ && normal_model_one != HOOKE) { - Emod[i][j] = Emod[j][i] = normal_coeffs_one[0]; - poiss[i][j] = poiss[j][i] = normal_coeffs_one[2]; - normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = - FOURTHIRDS*mix_stiffnessE(Emod[i][j],Emod[i][j], - poiss[i][j],poiss[i][j]); - } else { - normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = normal_coeffs_one[0]; - } - if ((normal_model_one == JKR) || (normal_model_one == DMT)) - normal_coeffs[i][j][3] = normal_coeffs[j][i][3] = normal_coeffs_one[3]; - - damping_model[i][j] = damping_model[j][i] = damping_model_one; - - tangential_model[i][j] = tangential_model[j][i] = tangential_model_one; - if (tangential_coeffs_one[0] == -1) { - tangential_coeffs[i][j][0] = tangential_coeffs[j][i][0] = - 8*mix_stiffnessG(Emod[i][j],Emod[i][j],poiss[i][j],poiss[i][j]); - } else { - tangential_coeffs[i][j][0] = tangential_coeffs[j][i][0] = - tangential_coeffs_one[0]; - } - for (int k = 1; k < 3; k++) - tangential_coeffs[i][j][k] = tangential_coeffs[j][i][k] = - tangential_coeffs_one[k]; - - roll_model[i][j] = roll_model[j][i] = roll_model_one; - if (roll_model_one != ROLL_NONE) - for (int k = 0; k < 3; k++) - roll_coeffs[i][j][k] = roll_coeffs[j][i][k] = roll_coeffs_one[k]; - - twist_model[i][j] = twist_model[j][i] = twist_model_one; - if (twist_model_one != TWIST_NONE && twist_model_one != TWIST_MARSHALL) - for (int k = 0; k < 3; k++) - twist_coeffs[i][j][k] = twist_coeffs[j][i][k] = twist_coeffs_one[k]; - - cutoff_type[i][j] = cutoff_type[j][i] = cutoff_one; - - setflag[i][j] = 1; - count++; - } - } - - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); -} - -/* ---------------------------------------------------------------------- - init specific to this pair style -------------------------------------------------------------------------- */ - -void PairGranular::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 = normal_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 = i; j <= atom->ntypes; j++) - if (normal_model[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; - } - } - for (int i = 1; i <= atom->ntypes; i++) - for (int j = i; j <= atom->ntypes; j++) - if (tangential_model[i][j] == TANGENTIAL_MINDLIN_RESCALE || - tangential_model[i][j] == TANGENTIAL_MINDLIN_RESCALE_FORCE) { - size_history += 1; - roll_history_index += 1; - twist_history_index += 1; - nondefault_history_transfer = 1; - history_transfer_factors = new int[size_history]; - for (int ii = 0; ii < size_history; ++ii) - history_transfer_factors[ii] = -1; - history_transfer_factors[3] = 1; - break; - } - - 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 and first init, create Fix to store history - // it replaces FixDummy, created in the constructor - // this is so its order in the fix list is preserved - - if (use_history && fix_history == nullptr) { - modify->replace_fix("NEIGH_HISTORY_GRANULAR_DUMMY","NEIGH_HISTORY_GRANULAR" - " all NEIGH_HISTORY " + std::to_string(size_history),1); - int ifix = modify->find_fix("NEIGH_HISTORY_GRANULAR"); - fix_history = (FixNeighHistory *) modify->fix[ifix]; - 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 = nullptr; - 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)); - onerad_dynamic[i] = radmax; - } - if (idep >= 0) { - itype = i; - double radmax = *((double *) modify->fix[idep]->extract("radius",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 (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_GRANULAR"); - 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 PairGranular::init_one(int i, int j) -{ - double cutoff=0.0; - - if (setflag[i][j] == 0) { - if ((normal_model[i][i] != normal_model[j][j]) || - (damping_model[i][i] != damping_model[j][j]) || - (tangential_model[i][i] != tangential_model[j][j]) || - (roll_model[i][i] != roll_model[j][j]) || - (twist_model[i][i] != twist_model[j][j])) { - - char str[512]; - sprintf(str,"Granular pair style functional forms are different, " - "cannot mix coefficients for types %d and %d. \n" - "This combination must be set explicitly " - "via pair_coeff command",i,j); - error->one(FLERR,str); - } - - if (normal_model[i][j] == HERTZ || normal_model[i][j] == HOOKE) - normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = - mix_geom(normal_coeffs[i][i][0], normal_coeffs[j][j][0]); - else - normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = - mix_stiffnessE(Emod[i][i], Emod[j][j], poiss[i][i], poiss[j][j]); - - normal_coeffs[i][j][1] = normal_coeffs[j][i][1] = - mix_geom(normal_coeffs[i][i][1], normal_coeffs[j][j][1]); - if ((normal_model[i][j] == JKR) || (normal_model[i][j] == 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] = tangential_coeffs[j][i][k] = - mix_geom(tangential_coeffs[i][i][k], tangential_coeffs[j][j][k]); - - if (roll_model[i][j] != 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_model[i][j] != TWIST_NONE && twist_model[i][j] != 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. - - double pulloff; - - if (cutoff_type[i][j] < 0 && cutoff_global < 0) { - if (((maxrad_dynamic[i] > 0.0) && (maxrad_dynamic[j] > 0.0)) || - ((maxrad_dynamic[i] > 0.0) && (maxrad_frozen[j] > 0.0)) || - // radius info about both i and j exist - ((maxrad_frozen[i] > 0.0) && (maxrad_dynamic[j] > 0.0))) { - cutoff = maxrad_dynamic[i]+maxrad_dynamic[j]; - pulloff = 0.0; - if (normal_model[i][j] == JKR) { - pulloff = pulloff_distance(maxrad_dynamic[i], maxrad_dynamic[j], i, j); - cutoff += pulloff; - } - - if (normal_model[i][j] == JKR) - pulloff = pulloff_distance(maxrad_frozen[i], maxrad_dynamic[j], i, j); - cutoff = MAX(cutoff, maxrad_frozen[i]+maxrad_dynamic[j]+pulloff); - - if (normal_model[i][j] == JKR) - pulloff = pulloff_distance(maxrad_dynamic[i], maxrad_frozen[j], i, j); - cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]+pulloff); - - } 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 if (cutoff_type[i][j] > 0) { - cutoff = cutoff_type[i][j]; - } else if (cutoff_global > 0) { - cutoff = cutoff_global; - } - - return cutoff; -} - -/* ---------------------------------------------------------------------- - proc 0 writes to restart file -------------------------------------------------------------------------- */ - -void PairGranular::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_model[i][j],sizeof(int),1,fp); - fwrite(&damping_model[i][j],sizeof(int),1,fp); - fwrite(&tangential_model[i][j],sizeof(int),1,fp); - fwrite(&roll_model[i][j],sizeof(int),1,fp); - fwrite(&twist_model[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(&cutoff_type[i][j],sizeof(double),1,fp); - } - } - } -} - -/* ---------------------------------------------------------------------- - proc 0 reads from restart file, bcasts -------------------------------------------------------------------------- */ - -void PairGranular::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) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error); - MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); - if (setflag[i][j]) { - if (me == 0) { - utils::sfread(FLERR,&normal_model[i][j],sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&damping_model[i][j],sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&tangential_model[i][j],sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&roll_model[i][j],sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&twist_model[i][j],sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,normal_coeffs[i][j],sizeof(double),4,fp,nullptr,error); - utils::sfread(FLERR,tangential_coeffs[i][j],sizeof(double),3,fp,nullptr,error); - utils::sfread(FLERR,roll_coeffs[i][j],sizeof(double),3,fp,nullptr,error); - utils::sfread(FLERR,twist_coeffs[i][j],sizeof(double),3,fp,nullptr,error); - utils::sfread(FLERR,&cutoff_type[i][j],sizeof(double),1,fp,nullptr,error); - } - MPI_Bcast(&normal_model[i][j],1,MPI_INT,0,world); - MPI_Bcast(&damping_model[i][j],1,MPI_INT,0,world); - MPI_Bcast(&tangential_model[i][j],1,MPI_INT,0,world); - MPI_Bcast(&roll_model[i][j],1,MPI_INT,0,world); - MPI_Bcast(&twist_model[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(&cutoff_type[i][j],1,MPI_DOUBLE,0,world); - } - } - } -} - -/* ---------------------------------------------------------------------- */ - -void PairGranular::reset_dt() -{ - dt = update->dt; -} - -/* ---------------------------------------------------------------------- */ - -double PairGranular::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,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; - double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3; - - double knfac, damp_normal, damp_normal_prefactor; - double k_tangential, damp_tangential; - double Fne, Ft, Fdamp, Fntot, Fncrit, 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; - - // rolling - double k_roll, damp_roll; - double rollmag; - double fr, fr1, fr2, fr3; - - // twisting - double k_twist, damp_twist, mu_twist; - double signtwist, magtwist, magtortwist, Mtcrit; - - double shrmag; - int jnum; - int *jlist; - double *history,*allhistory; - - int nall = atom->nlocal + atom->nghost; - if ((i >= nall) || (j >= nall)) - error->all(FLERR,"Not enough atoms for pair granular single function"); - - double *radius = atom->radius; - radi = radius[i]; - radj = radius[j]; - radsum = radi + radj; - Reff = (radsum > 0.0) ? radi*radj/radsum : 0.0; - - bool touchflag; - E = normal_coeffs[itype][jtype][0]; - if (normal_model[itype][jtype] == JKR) { - E *= THREEQUARTERS; - R2 = Reff*Reff; - coh = normal_coeffs[itype][jtype][3]; - a = cbrt(9.0*MY_PI*coh*R2/(4*E)); - delta_pulloff = a*a/Reff - 2*sqrt(MY_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; - - // 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 *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; - - delta = radsum - r; - dR = delta*Reff; - if (normal_model[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; - // in case sqrt(0) < 0 due to precision issues - sqrt1 = MAX(0, t0*(t1+2*t2)); - t4 = cbrt(t1+t2+THREEROOT3*MY_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*MY_PI*R2/(E*t6)); - a = INVROOT6*(t6 + sqrt(sqrt3)); - a2 = a*a; - knfac = normal_coeffs[itype][jtype][0]*a; - Fne = knfac*a2/Reff - MY_2PI*a2*sqrt(4*coh*E/(MY_PI*a)); - } else { - knfac = E; - Fne = knfac*delta; - a = sqrt(dR); - if (normal_model[itype][jtype] != HOOKE) { - Fne *= a; - knfac *= a; - } - if (normal_model[itype][jtype] == DMT) - Fne -= 4*MY_PI*normal_coeffs[itype][jtype][3]*Reff; - } - - if (damping_model[itype][jtype] == VELOCITY) { - damp_normal = 1; - } else if (damping_model[itype][jtype] == MASS_VELOCITY) { - damp_normal = meff; - } else if (damping_model[itype][jtype] == VISCOELASTIC) { - damp_normal = a*meff; - } else if (damping_model[itype][jtype] == TSUJI) { - damp_normal = sqrt(meff*knfac); - } else damp_normal = 0.0; - - damp_normal_prefactor = normal_coeffs[itype][jtype][1]*damp_normal; - Fdamp = -damp_normal_prefactor*vnnr; - - Fntot = Fne + Fdamp; - - jnum = list->numneigh[i]; - jlist = list->firstneigh[i]; - - if (use_history) { - if ((fix_history == nullptr) || (fix_history->firstvalue == nullptr)) - error->one(FLERR,"Pair granular single computation needs 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 - //**************************************** - - // For linear, mindlin, mindlin_rescale: - // history = cumulative tangential displacement - // - // For mindlin/force, mindlin_rescale/force: - // history = cumulative tangential elastic force - - // 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_model[itype][jtype] == JKR) { - F_pulloff = 3*MY_PI*coh*Reff; - Fncrit = fabs(Fne + 2*F_pulloff); - } else if (normal_model[itype][jtype] == DMT) { - F_pulloff = 4*MY_PI*coh*Reff; - Fncrit = fabs(Fne + 2*F_pulloff); - } else { - Fncrit = fabs(Fntot); - } - Fscrit = tangential_coeffs[itype][jtype][2] * Fncrit; - - //------------------------------ - // tangential forces - //------------------------------ - k_tangential = tangential_coeffs[itype][jtype][0]; - damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal_prefactor; - - if (tangential_history) { - if (tangential_model[itype][jtype] != TANGENTIAL_HISTORY) { - k_tangential *= a; - } - - shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + - history[2]*history[2]); - - // tangential forces = history + tangential velocity damping - if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY || - tangential_model[itype][jtype] == TANGENTIAL_MINDLIN || - tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) { - fs1 = -k_tangential*history[0] - damp_tangential*vtr1; - fs2 = -k_tangential*history[1] - damp_tangential*vtr2; - fs3 = -k_tangential*history[2] - damp_tangential*vtr3; - } else { - fs1 = history[0] - damp_tangential*vtr1; - fs2 = history[1] - damp_tangential*vtr2; - fs3 = history[2] - damp_tangential*vtr3; - } - - // rescale frictional forces if needed - fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); - 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; - } - - // classic pair gran/hooke (no history) - } else { - fs = damp_tangential*vrel; - if (vrel != 0.0) Ft = MIN(Fscrit,fs) / vrel; - else Ft = 0.0; - fs1 = -Ft*vtr1; - fs2 = -Ft*vtr2; - fs3 = -Ft*vtr3; - fs = Ft*vrel; - } - - //**************************************** - // rolling resistance - //**************************************** - - if ((roll_model[itype][jtype] != ROLL_NONE) - || (twist_model[itype][jtype] != TWIST_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; - - 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]); - - 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] * Fncrit; - - fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3); - if (fr > Frcrit) { - if (rollmag != 0.0) { - fr1 *= Frcrit/fr; - fr2 *= Frcrit/fr; - fr3 *= Frcrit/fr; - fr *= Frcrit/fr; - } else fr1 = fr2 = fr3 = fr = 0.0; - } - } else fr1 = fr2 = fr3 = fr = 0.0; - - //**************************************** - // twisting torque, including history effects - //**************************************** - - if (twist_model[itype][jtype] != TWIST_NONE) { - // omega_T (eq 29 of Marshall) - magtwist = relrot1*nx + relrot2*ny + relrot3*nz; - if (twist_model[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*tangential_coeffs[itype][jtype][2];; - } else { - k_twist = twist_coeffs[itype][jtype][0]; - damp_twist = twist_coeffs[itype][jtype][1]; - mu_twist = twist_coeffs[itype][jtype][2]; - } - // M_t torque (eq 30) - magtortwist = -k_twist*history[twist_history_index] - damp_twist*magtwist; - signtwist = (magtwist > 0) - (magtwist < 0); - Mtcrit = mu_twist*Fncrit; // critical torque (eq 44) - if (fabs(magtortwist) > Mtcrit) { - magtortwist = -Mtcrit * signtwist; // eq 34 - } else magtortwist = 0.0; - } else magtortwist = 0.0; - - // set force and return no energy - - fforce = Fntot*rinv; - - // 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; - svector[9] = delx; - svector[10] = dely; - svector[11] = delz; - return 0.0; -} - -/* ---------------------------------------------------------------------- */ - -int PairGranular::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 PairGranular::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 PairGranular::memory_usage() -{ - double bytes = (double)nmax * sizeof(double); - return bytes; -} - -/* ---------------------------------------------------------------------- - mixing of Young's modulus (E) -------------------------------------------------------------------------- */ - -double PairGranular::mix_stiffnessE(double Eii, double Ejj, - double poisii, double poisjj) -{ - return 1/((1-poisii*poisii)/Eii+(1-poisjj*poisjj)/Ejj); -} - -/* ---------------------------------------------------------------------- - mixing of shear modulus (G) ------------------------------------------------------------------------- */ - -double PairGranular::mix_stiffnessG(double Eii, double Ejj, - double poisii, double poisjj) -{ - return 1/((2*(2-poisii)*(1+poisii)/Eii) + (2*(2-poisjj)*(1+poisjj)/Ejj)); -} - -/* ---------------------------------------------------------------------- - mixing of everything else -------------------------------------------------------------------------- */ - -double PairGranular::mix_geom(double valii, double valjj) -{ - return sqrt(valii*valjj); -} - - -/* ---------------------------------------------------------------------- - compute pull-off distance (beyond contact) for a given radius and atom type -------------------------------------------------------------------------- */ - -double PairGranular::pulloff_distance(double radi, double radj, - int itype, int jtype) -{ - double E, coh, a, Reff; - Reff = radi*radj/(radi+radj); - if (Reff <= 0) return 0; - coh = normal_coeffs[itype][jtype][3]; - E = normal_coeffs[itype][jtype][0]*THREEQUARTERS; - a = cbrt(9*MY_PI*coh*Reff*Reff/(4*E)); - return a*a/Reff - 2*sqrt(MY_PI*coh*a/E); -} - -/* ---------------------------------------------------------------------- - transfer history during fix/neigh/history exchange - only needed if any history entries i-j are not just negative of j-i entries -------------------------------------------------------------------------- */ - -void PairGranular::transfer_history(double* source, double* target) -{ - for (int i = 0; i < size_history; i++) - target[i] = history_transfer_factors[i]*source[i]; -} From 6bdf0138acbeb53101bc7ac3b2e42cfd995470da Mon Sep 17 00:00:00 2001 From: Joel Clemmer Date: Mon, 5 Apr 2021 11:50:53 -0600 Subject: [PATCH 08/15] Typos in documentation --- doc/src/pair_gran.rst | 2 +- doc/src/pair_granular.rst | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/src/pair_gran.rst b/doc/src/pair_gran.rst index 7d6189fda1..5bccdfd8b4 100644 --- a/doc/src/pair_gran.rst +++ b/doc/src/pair_gran.rst @@ -219,7 +219,7 @@ potential. If two particles are moving away from each other while in contact, there is a possibility that the particles could experience an effective attractive -force due to damping. If the *limit_damping* keyword is used, this fix +force due to damping. If the *limit_damping* keyword is used, this option will zero out the normal component of the force if there is an effective attractive force. diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index e0a33f31d0..a9ca437370 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -625,7 +625,7 @@ Finally, the twisting torque on each particle is given by: If two particles are moving away from each other while in contact, there is a possibility that the particles could experience an effective attractive -force due to damping. If the optional *limit_damping* keyword is used, this fix +force due to damping. If the optional *limit_damping* keyword is used, this option will zero out the normal component of the force if there is an effective attractive force. This keyword cannot be used with the JKR or DMT models. @@ -667,7 +667,7 @@ atom types. If two particles are moving away from each other while in contact, there is a possibility that the particles could experience an effective attractive -force due to damping. If the *no_attraction* keyword is used, this fix +force due to damping. If the *limit_damping* keyword is used, this option will zero out the normal component of the force if there is an effective attractive force. This keyword cannot be used with the JKR or DMT models. From f323fb29b3cc6e96a643b9de94c73472c3f9bccd Mon Sep 17 00:00:00 2001 From: Joel Clemmer Date: Mon, 5 Apr 2021 11:53:49 -0600 Subject: [PATCH 09/15] Reverting no_attraction option in wall gran region docu --- doc/src/fix_wall_gran_region.rst | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/doc/src/fix_wall_gran_region.rst b/doc/src/fix_wall_gran_region.rst index b7a9ff6bd7..35fe1fab72 100644 --- a/doc/src/fix_wall_gran_region.rst +++ b/doc/src/fix_wall_gran_region.rst @@ -36,14 +36,12 @@ Syntax * wallstyle = region (see :doc:`fix wall/gran ` for options for other kinds of walls) * region-ID = region whose boundary will act as wall -* keyword = *contacts* or *no_attraction* +* keyword = *contacts* .. parsed-literal:: *contacts* value = none generate contact information for each particle - *no_attraction* value = none - turn off possibility of attractive interactions Examples """""""" From 80a65150c22c186e230cb7d2902f2951065dc83a Mon Sep 17 00:00:00 2001 From: Joel Clemmer Date: Tue, 6 Apr 2021 13:08:48 -0600 Subject: [PATCH 10/15] Patching uninitialized values, replacing lines in documentation --- doc/src/fix_wall_gran.rst | 2 + src/GRANULAR/fix_wall_gran.cpp | 76 ++++++++++++------------ src/GRANULAR/pair_gran_hertz_history.cpp | 8 ++- src/GRANULAR/pair_granular.cpp | 2 +- 4 files changed, 49 insertions(+), 39 deletions(-) diff --git a/doc/src/fix_wall_gran.rst b/doc/src/fix_wall_gran.rst index 02d52ef1a6..95a4b5d818 100644 --- a/doc/src/fix_wall_gran.rst +++ b/doc/src/fix_wall_gran.rst @@ -177,6 +177,8 @@ the clockwise direction for *vshear* > 0 or counter-clockwise for *vshear* < 0. In this case, *vshear* is the tangential velocity of the wall at whatever *radius* has been defined. + +Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" This fix writes the shear friction state of atoms interacting with the diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index d746923ab0..f9840c9d5e 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -71,6 +71,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Fix wall/gran requires atom style sphere"); create_attribute = 1; + limit_damping = 0; // set interaction style // disable bonded/history option for now @@ -91,7 +92,6 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : // wall/particle coefficients int iarg; - if (pairstyle != GRANULAR) { size_history = 3; if (narg < 11) error->all(FLERR,"Illegal fix wall/gran command"); @@ -322,14 +322,14 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : if (normal_model == JKR) size_history += 1; if (tangential_model == TANGENTIAL_MINDLIN_RESCALE || tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) size_history += 1; - } - if (limit_damping and normal_model == JKR) - error->all(FLERR,"Illegal pair_coeff command, " - "cannot limit damping with JRK model"); - if (limit_damping and normal_model == DMT) - error->all(FLERR,"Illegal pair_coeff command, " - "Cannot limit damping with DMT model"); + if (limit_damping and normal_model == JKR) + error->all(FLERR,"Illegal pair_coeff command, " + "cannot limit damping with JRK model"); + if (limit_damping and normal_model == DMT) + error->all(FLERR,"Illegal pair_coeff command, " + "Cannot limit damping with DMT model"); + } // wallstyle args @@ -504,37 +504,39 @@ void FixWallGran::init() if (modify->fix[i]->rigid_flag) break; if (i < modify->nfix) fix_rigid = modify->fix[i]; - tangential_history_index = 0; - 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; + if(pairstyle == GRANULAR) { + tangential_history_index = 0; + if (roll_history) { + if (tangential_history) roll_history_index = 3; + else roll_history_index = 0; } - else{ - if (roll_history) twist_history_index = 3; - else twist_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; + } + } + if (normal_model == JKR) { + tangential_history_index += 1; + roll_history_index += 1; + twist_history_index += 1; + } + if (tangential_model == TANGENTIAL_MINDLIN_RESCALE || + tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) { + roll_history_index += 1; + twist_history_index += 1; + } + + if (damping_model == TSUJI) { + double cor = normal_coeffs[1]; + normal_coeffs[1] = 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); } - } - if (normal_model == JKR) { - tangential_history_index += 1; - roll_history_index += 1; - twist_history_index += 1; - } - if (tangential_model == TANGENTIAL_MINDLIN_RESCALE || - tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) { - roll_history_index += 1; - twist_history_index += 1; - } - - if (damping_model == TSUJI) { - double cor = normal_coeffs[1]; - normal_coeffs[1] = 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); } } diff --git a/src/GRANULAR/pair_gran_hertz_history.cpp b/src/GRANULAR/pair_gran_hertz_history.cpp index e21ec727d6..be2a3c2558 100644 --- a/src/GRANULAR/pair_gran_hertz_history.cpp +++ b/src/GRANULAR/pair_gran_hertz_history.cpp @@ -278,7 +278,7 @@ void PairGranHertzHistory::compute(int eflag, int vflag) void PairGranHertzHistory::settings(int narg, char **arg) { - if (narg != 6) error->all(FLERR,"Illegal pair_style command"); + if (narg != 6 and narg != 7) error->all(FLERR,"Illegal pair_style command"); kn = utils::numeric(FLERR,arg[0],false,lmp); if (strcmp(arg[1],"NULL") == 0) kt = kn * 2.0/7.0; @@ -296,6 +296,12 @@ void PairGranHertzHistory::settings(int narg, char **arg) xmu < 0.0 || xmu > 10000.0 || dampflag < 0 || dampflag > 1) error->all(FLERR,"Illegal pair_style command"); + limit_damping = 0; + if (narg == 7) { + if (strcmp(arg[6], "limit_damping") == 0) limit_damping = 1; + else error->all(FLERR,"Illegal pair_style command"); + } + // convert Kn and Kt from pressure units to force/distance^2 kn /= force->nktv2p; diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index ac98a6e27c..653d5c70a1 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -1038,7 +1038,7 @@ void PairGranular::coeff(int narg, char **arg) cutoff_type[i][j] = cutoff_type[j][i] = cutoff_one; - limit_damping[i][j] = ld_flag; + limit_damping[i][j] = limit_damping[j][i] = ld_flag; setflag[i][j] = 1; count++; From 8e43d58faba0bcb520bdaa0e21c9a5f1eb601017 Mon Sep 17 00:00:00 2001 From: Joel Clemmer Date: Tue, 6 Apr 2021 13:48:56 -0600 Subject: [PATCH 11/15] Switching logical operators to match preferred style --- src/GRANULAR/fix_wall_gran.cpp | 12 ++++++------ src/GRANULAR/pair_gran_hertz_history.cpp | 6 +++--- src/GRANULAR/pair_gran_hooke.cpp | 4 ++-- src/GRANULAR/pair_gran_hooke_history.cpp | 6 +++--- src/GRANULAR/pair_granular.cpp | 8 ++++---- src/KOKKOS/pair_gran_hooke_history_kokkos.cpp | 2 +- 6 files changed, 19 insertions(+), 19 deletions(-) diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index f9840c9d5e..11cf0a9bc9 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -323,10 +323,10 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : if (tangential_model == TANGENTIAL_MINDLIN_RESCALE || tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) size_history += 1; - if (limit_damping and normal_model == JKR) + if (limit_damping && normal_model == JKR) error->all(FLERR,"Illegal pair_coeff command, " "cannot limit damping with JRK model"); - if (limit_damping and normal_model == DMT) + if (limit_damping && normal_model == DMT) error->all(FLERR,"Illegal pair_coeff command, " "Cannot limit damping with DMT model"); } @@ -794,7 +794,7 @@ void FixWallGran::hooke(double rsq, double dx, double dy, double dz, damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radius-r)*rinv - damp; - if (limit_damping and ccel < 0.0) ccel = 0.0; + if (limit_damping && ccel < 0.0) ccel = 0.0; // relative velocities @@ -887,7 +887,7 @@ void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz, damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radius-r)*rinv - damp; - if (limit_damping and ccel < 0.0) ccel = 0.0; + if (limit_damping && ccel < 0.0) ccel = 0.0; // relative velocities @@ -1019,7 +1019,7 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz, if (rwall == 0.0) polyhertz = sqrt((radius-r)*radius); else polyhertz = sqrt((radius-r)*radius*rwall/(rwall+radius)); ccel *= polyhertz; - if (limit_damping and ccel < 0.0) ccel = 0.0; + if (limit_damping && ccel < 0.0) ccel = 0.0; // relative velocities @@ -1214,7 +1214,7 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz, Fdamp = -damp_normal_prefactor*vnnr; Fntot = Fne + Fdamp; - if (limit_damping and Fntot < 0.0) Fntot = 0.0; + if (limit_damping && Fntot < 0.0) Fntot = 0.0; //**************************************** // tangential force, including history effects diff --git a/src/GRANULAR/pair_gran_hertz_history.cpp b/src/GRANULAR/pair_gran_hertz_history.cpp index be2a3c2558..e40a9be2bd 100644 --- a/src/GRANULAR/pair_gran_hertz_history.cpp +++ b/src/GRANULAR/pair_gran_hertz_history.cpp @@ -181,7 +181,7 @@ void PairGranHertzHistory::compute(int eflag, int vflag) ccel = kn*(radsum-r)*rinv - damp; polyhertz = sqrt((radsum-r)*radi*radj / radsum); ccel *= polyhertz; - if (limit_damping and ccel < 0.0) ccel = 0.0; + if (limit_damping && ccel < 0.0) ccel = 0.0; // relative velocities @@ -278,7 +278,7 @@ void PairGranHertzHistory::compute(int eflag, int vflag) void PairGranHertzHistory::settings(int narg, char **arg) { - if (narg != 6 and narg != 7) error->all(FLERR,"Illegal pair_style command"); + if (narg != 6 && narg != 7) error->all(FLERR,"Illegal pair_style command"); kn = utils::numeric(FLERR,arg[0],false,lmp); if (strcmp(arg[1],"NULL") == 0) kt = kn * 2.0/7.0; @@ -395,7 +395,7 @@ double PairGranHertzHistory::single(int i, int j, int /*itype*/, int /*jtype*/, ccel = kn*(radsum-r)*rinv - damp; polyhertz = sqrt((radsum-r)*radi*radj / radsum); ccel *= polyhertz; - if (limit_damping and ccel < 0.0) ccel = 0.0; + if (limit_damping && ccel < 0.0) ccel = 0.0; // relative velocities diff --git a/src/GRANULAR/pair_gran_hooke.cpp b/src/GRANULAR/pair_gran_hooke.cpp index c28bbd007f..0f45fadca9 100644 --- a/src/GRANULAR/pair_gran_hooke.cpp +++ b/src/GRANULAR/pair_gran_hooke.cpp @@ -158,7 +158,7 @@ void PairGranHooke::compute(int eflag, int vflag) damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; - if (limit_damping and ccel < 0.0) ccel = 0.0; + if (limit_damping && ccel < 0.0) ccel = 0.0; // relative velocities @@ -300,7 +300,7 @@ double PairGranHooke::single(int i, int j, int /*itype*/, int /*jtype*/, double damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; - if (limit_damping and ccel < 0.0) ccel = 0.0; + if (limit_damping && ccel < 0.0) ccel = 0.0; // relative velocities diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp index 4f7a6e20ec..48c30c761e 100644 --- a/src/GRANULAR/pair_gran_hooke_history.cpp +++ b/src/GRANULAR/pair_gran_hooke_history.cpp @@ -237,7 +237,7 @@ void PairGranHookeHistory::compute(int eflag, int vflag) damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; - if (limit_damping and ccel < 0.0) ccel = 0.0; + if (limit_damping && ccel < 0.0) ccel = 0.0; // relative velocities @@ -357,7 +357,7 @@ void PairGranHookeHistory::allocate() void PairGranHookeHistory::settings(int narg, char **arg) { - if (narg != 6 and narg != 7) error->all(FLERR,"Illegal pair_style command"); + if (narg != 6 && narg != 7) error->all(FLERR,"Illegal pair_style command"); kn = utils::numeric(FLERR,arg[0],false,lmp); if (strcmp(arg[1],"NULL") == 0) kt = kn * 2.0/7.0; @@ -693,7 +693,7 @@ double PairGranHookeHistory::single(int i, int j, int /*itype*/, int /*jtype*/, damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; - if(limit_damping and ccel < 0.0) ccel = 0.0; + if(limit_damping && ccel < 0.0) ccel = 0.0; // relative velocities diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 653d5c70a1..e787e305ff 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -371,7 +371,7 @@ void PairGranular::compute(int eflag, int vflag) Fdamp = -damp_normal_prefactor*vnnr; Fntot = Fne + Fdamp; - if (limit_damping[itype][jtype] and Fntot < 0.0) Fntot = 0.0; + if (limit_damping[itype][jtype] && Fntot < 0.0) Fntot = 0.0; //**************************************** // tangential force, including history effects @@ -988,11 +988,11 @@ void PairGranular::coeff(int narg, char **arg) 27.467*powint(cor,4)-18.022*powint(cor,5)+4.8218*powint(cor,6); } else damp = normal_coeffs_one[1]; - if (ld_flag and normal_model_one == JKR) + if (ld_flag && normal_model_one == JKR) error->all(FLERR,"Illegal pair_coeff command, " "Cannot limit damping with JKR model"); - if (ld_flag and normal_model_one == DMT) + if (ld_flag && normal_model_one == DMT) error->all(FLERR,"Illegal pair_coeff command, " "Cannot limit damping with DMT model"); @@ -1547,7 +1547,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, Fdamp = -damp_normal_prefactor*vnnr; Fntot = Fne + Fdamp; - if (limit_damping[itype][jtype] and Fntot < 0.0) Fntot = 0.0; + if (limit_damping[itype][jtype] && Fntot < 0.0) Fntot = 0.0; jnum = list->numneigh[i]; jlist = list->firstneigh[i]; diff --git a/src/KOKKOS/pair_gran_hooke_history_kokkos.cpp b/src/KOKKOS/pair_gran_hooke_history_kokkos.cpp index 8d853b7dae..95be9843b9 100644 --- a/src/KOKKOS/pair_gran_hooke_history_kokkos.cpp +++ b/src/KOKKOS/pair_gran_hooke_history_kokkos.cpp @@ -392,7 +392,7 @@ void PairGranHookeHistoryKokkos::operator()(TagPairGranHookeHistoryC F_FLOAT damp = meff*gamman*vnnr*rsqinv; F_FLOAT ccel = kn*(radsum-r)*rinv - damp; - if(limit_damping & ccel < 0.0) ccel = 0.0; + if(limit_damping && ccel < 0.0) ccel = 0.0; // relative velocities From c44b8f18ee49326cf63bf14410bb1c4493cfb700 Mon Sep 17 00:00:00 2001 From: Joel Clemmer Date: Tue, 6 Apr 2021 13:54:11 -0600 Subject: [PATCH 12/15] Adding explicit parentheses to logical operations --- src/GRANULAR/fix_wall_gran.cpp | 8 ++++---- src/GRANULAR/pair_gran_hertz_history.cpp | 4 ++-- src/GRANULAR/pair_gran_hooke.cpp | 4 ++-- src/GRANULAR/pair_gran_hooke_history.cpp | 4 ++-- src/GRANULAR/pair_granular.cpp | 4 ++-- src/KOKKOS/pair_gran_hooke_history_kokkos.cpp | 2 +- 6 files changed, 13 insertions(+), 13 deletions(-) diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index 11cf0a9bc9..0611a54efd 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -794,7 +794,7 @@ void FixWallGran::hooke(double rsq, double dx, double dy, double dz, damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radius-r)*rinv - damp; - if (limit_damping && ccel < 0.0) ccel = 0.0; + if (limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities @@ -887,7 +887,7 @@ void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz, damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radius-r)*rinv - damp; - if (limit_damping && ccel < 0.0) ccel = 0.0; + if (limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities @@ -1019,7 +1019,7 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz, if (rwall == 0.0) polyhertz = sqrt((radius-r)*radius); else polyhertz = sqrt((radius-r)*radius*rwall/(rwall+radius)); ccel *= polyhertz; - if (limit_damping && ccel < 0.0) ccel = 0.0; + if (limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities @@ -1214,7 +1214,7 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz, Fdamp = -damp_normal_prefactor*vnnr; Fntot = Fne + Fdamp; - if (limit_damping && Fntot < 0.0) Fntot = 0.0; + if (limit_damping && (Fntot < 0.0)) Fntot = 0.0; //**************************************** // tangential force, including history effects diff --git a/src/GRANULAR/pair_gran_hertz_history.cpp b/src/GRANULAR/pair_gran_hertz_history.cpp index e40a9be2bd..56e9691fad 100644 --- a/src/GRANULAR/pair_gran_hertz_history.cpp +++ b/src/GRANULAR/pair_gran_hertz_history.cpp @@ -181,7 +181,7 @@ void PairGranHertzHistory::compute(int eflag, int vflag) ccel = kn*(radsum-r)*rinv - damp; polyhertz = sqrt((radsum-r)*radi*radj / radsum); ccel *= polyhertz; - if (limit_damping && ccel < 0.0) ccel = 0.0; + if (limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities @@ -395,7 +395,7 @@ double PairGranHertzHistory::single(int i, int j, int /*itype*/, int /*jtype*/, ccel = kn*(radsum-r)*rinv - damp; polyhertz = sqrt((radsum-r)*radi*radj / radsum); ccel *= polyhertz; - if (limit_damping && ccel < 0.0) ccel = 0.0; + if (limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities diff --git a/src/GRANULAR/pair_gran_hooke.cpp b/src/GRANULAR/pair_gran_hooke.cpp index 0f45fadca9..eb5ccec680 100644 --- a/src/GRANULAR/pair_gran_hooke.cpp +++ b/src/GRANULAR/pair_gran_hooke.cpp @@ -158,7 +158,7 @@ void PairGranHooke::compute(int eflag, int vflag) damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; - if (limit_damping && ccel < 0.0) ccel = 0.0; + if (limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities @@ -300,7 +300,7 @@ double PairGranHooke::single(int i, int j, int /*itype*/, int /*jtype*/, double damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; - if (limit_damping && ccel < 0.0) ccel = 0.0; + if (limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp index 48c30c761e..0f9682a133 100644 --- a/src/GRANULAR/pair_gran_hooke_history.cpp +++ b/src/GRANULAR/pair_gran_hooke_history.cpp @@ -237,7 +237,7 @@ void PairGranHookeHistory::compute(int eflag, int vflag) damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; - if (limit_damping && ccel < 0.0) ccel = 0.0; + if (limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities @@ -693,7 +693,7 @@ double PairGranHookeHistory::single(int i, int j, int /*itype*/, int /*jtype*/, damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; - if(limit_damping && ccel < 0.0) ccel = 0.0; + if(limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index e787e305ff..dc367dcc0a 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -371,7 +371,7 @@ void PairGranular::compute(int eflag, int vflag) Fdamp = -damp_normal_prefactor*vnnr; Fntot = Fne + Fdamp; - if (limit_damping[itype][jtype] && Fntot < 0.0) Fntot = 0.0; + if (limit_damping[itype][jtype] && (Fntot < 0.0)) Fntot = 0.0; //**************************************** // tangential force, including history effects @@ -1547,7 +1547,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, Fdamp = -damp_normal_prefactor*vnnr; Fntot = Fne + Fdamp; - if (limit_damping[itype][jtype] && Fntot < 0.0) Fntot = 0.0; + if (limit_damping[itype][jtype] && (Fntot < 0.0)) Fntot = 0.0; jnum = list->numneigh[i]; jlist = list->firstneigh[i]; diff --git a/src/KOKKOS/pair_gran_hooke_history_kokkos.cpp b/src/KOKKOS/pair_gran_hooke_history_kokkos.cpp index 95be9843b9..01ac1ea54f 100644 --- a/src/KOKKOS/pair_gran_hooke_history_kokkos.cpp +++ b/src/KOKKOS/pair_gran_hooke_history_kokkos.cpp @@ -392,7 +392,7 @@ void PairGranHookeHistoryKokkos::operator()(TagPairGranHookeHistoryC F_FLOAT damp = meff*gamman*vnnr*rsqinv; F_FLOAT ccel = kn*(radsum-r)*rinv - damp; - if(limit_damping && ccel < 0.0) ccel = 0.0; + if(limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities From da56b9de569ac6546df64872283ae6f2eadde983 Mon Sep 17 00:00:00 2001 From: Joel Clemmer Date: Wed, 7 Apr 2021 07:48:43 -0600 Subject: [PATCH 13/15] Adding user-omp support --- src/USER-OMP/pair_gran_hertz_history_omp.cpp | 1 + src/USER-OMP/pair_gran_hooke_history_omp.cpp | 1 + src/USER-OMP/pair_gran_hooke_omp.cpp | 1 + 3 files changed, 3 insertions(+) diff --git a/src/USER-OMP/pair_gran_hertz_history_omp.cpp b/src/USER-OMP/pair_gran_hertz_history_omp.cpp index 2187dcbe90..200075a316 100644 --- a/src/USER-OMP/pair_gran_hertz_history_omp.cpp +++ b/src/USER-OMP/pair_gran_hertz_history_omp.cpp @@ -224,6 +224,7 @@ void PairGranHertzHistoryOMP::eval(int iifrom, int iito, ThrData * const thr) ccel = kn*(radsum-r)*rinv - damp; polyhertz = sqrt((radsum-r)*radi*radj / radsum); ccel *= polyhertz; + if (limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities diff --git a/src/USER-OMP/pair_gran_hooke_history_omp.cpp b/src/USER-OMP/pair_gran_hooke_history_omp.cpp index 7894d1e963..0ba9607bd9 100644 --- a/src/USER-OMP/pair_gran_hooke_history_omp.cpp +++ b/src/USER-OMP/pair_gran_hooke_history_omp.cpp @@ -223,6 +223,7 @@ void PairGranHookeHistoryOMP::eval(int iifrom, int iito, ThrData * const thr) damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; + if (limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities diff --git a/src/USER-OMP/pair_gran_hooke_omp.cpp b/src/USER-OMP/pair_gran_hooke_omp.cpp index 08ae07c6f5..451b074fa9 100644 --- a/src/USER-OMP/pair_gran_hooke_omp.cpp +++ b/src/USER-OMP/pair_gran_hooke_omp.cpp @@ -197,6 +197,7 @@ void PairGranHookeOMP::eval(int iifrom, int iito, ThrData * const thr) damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; + if (limit_damping && (ccel < 0.0)) ccel = 0.0; // relative velocities From f072289ac1ad17fca3aaf23046ca8ca0aa1bb74c Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 7 Apr 2021 14:52:36 -0400 Subject: [PATCH 14/15] need to initialize limit_damping array to NULL --- src/GRANULAR/pair_granular.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index dc367dcc0a..faeaaacd88 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -80,6 +80,8 @@ PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp) onerad_frozen = nullptr; maxrad_dynamic = nullptr; maxrad_frozen = nullptr; + + limit_damping = nullptr; history_transfer_factors = nullptr; From ea8277ce87d8e5b7fe5be83117f9e4025687662f Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 7 Apr 2021 14:59:33 -0400 Subject: [PATCH 15/15] whitespace fixes --- doc/src/fix_wall_gran.rst | 4 ++-- doc/src/fix_wall_gran_region.rst | 4 ++-- doc/src/pair_gran.rst | 8 ++++---- doc/src/pair_granular.rst | 12 ++++++------ src/GRANULAR/fix_wall_gran.cpp | 14 +++++++------- src/GRANULAR/pair_gran_hooke_history.cpp | 2 +- src/GRANULAR/pair_granular.cpp | 10 +++++----- 7 files changed, 27 insertions(+), 27 deletions(-) diff --git a/doc/src/fix_wall_gran.rst b/doc/src/fix_wall_gran.rst index 95a4b5d818..57bec679e4 100644 --- a/doc/src/fix_wall_gran.rst +++ b/doc/src/fix_wall_gran.rst @@ -96,8 +96,8 @@ Specifically, delta = radius - r = overlap of particle with wall, m_eff = mass of particle, and the effective radius of contact = RiRj/Ri+Rj is set to the radius of the particle. -The parameters *Kn*\ , *Kt*\ , *gamma_n*, *gamma_t*, *xmu*, *dampflag*, -and the optional keyword *limit_damping* +The parameters *Kn*\ , *Kt*\ , *gamma_n*, *gamma_t*, *xmu*, *dampflag*, +and the optional keyword *limit_damping* have the same meaning and units as those specified with the :doc:`pair_style gran/\* ` commands. This means a NULL can be used for either *Kt* or *gamma_t* as described on that page. If a diff --git a/doc/src/fix_wall_gran_region.rst b/doc/src/fix_wall_gran_region.rst index 35fe1fab72..0a252b162a 100644 --- a/doc/src/fix_wall_gran_region.rst +++ b/doc/src/fix_wall_gran_region.rst @@ -181,8 +181,8 @@ radius - r = overlap of particle with wall, m_eff = mass of particle, and the effective radius of contact is just the radius of the particle. -The parameters *Kn*\ , *Kt*\ , *gamma_n*, *gamma_t*, *xmu*, *dampflag*, -and the optional keyword *limit_damping* +The parameters *Kn*\ , *Kt*\ , *gamma_n*, *gamma_t*, *xmu*, *dampflag*, +and the optional keyword *limit_damping* have the same meaning and units as those specified with the :doc:`pair_style gran/\* ` commands. This means a NULL can be used for either *Kt* or *gamma_t* as described on that page. If a diff --git a/doc/src/pair_gran.rst b/doc/src/pair_gran.rst index 5bccdfd8b4..fbcacb5c76 100644 --- a/doc/src/pair_gran.rst +++ b/doc/src/pair_gran.rst @@ -217,11 +217,11 @@ potential is used as a sub-style of :doc:`pair_style hybrid `, then pair_coeff command to determine which atoms interact via a granular potential. -If two particles are moving away from each other while in contact, there -is a possibility that the particles could experience an effective attractive -force due to damping. If the *limit_damping* keyword is used, this option +If two particles are moving away from each other while in contact, there +is a possibility that the particles could experience an effective attractive +force due to damping. If the *limit_damping* keyword is used, this option will zero out the normal component of the force if there is an effective -attractive force. +attractive force. ---------- diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index a9ca437370..432fd29584 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -623,9 +623,9 @@ Finally, the twisting torque on each particle is given by: ---------- -If two particles are moving away from each other while in contact, there -is a possibility that the particles could experience an effective attractive -force due to damping. If the optional *limit_damping* keyword is used, this option +If two particles are moving away from each other while in contact, there +is a possibility that the particles could experience an effective attractive +force due to damping. If the optional *limit_damping* keyword is used, this option will zero out the normal component of the force if there is an effective attractive force. This keyword cannot be used with the JKR or DMT models. @@ -665,9 +665,9 @@ then LAMMPS will use that cutoff for the specified atom type combination, and automatically set pairwise cutoffs for the remaining atom types. -If two particles are moving away from each other while in contact, there -is a possibility that the particles could experience an effective attractive -force due to damping. If the *limit_damping* keyword is used, this option +If two particles are moving away from each other while in contact, there +is a possibility that the particles could experience an effective attractive +force due to damping. If the *limit_damping* keyword is used, this option will zero out the normal component of the force if there is an effective attractive force. This keyword cannot be used with the JKR or DMT models. diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index 0611a54efd..51dca15e7b 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -119,12 +119,12 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : kt /= force->nktv2p; } iarg = 10; - + if (strcmp(arg[iarg],"limit_damping") == 0) { limit_damping = 1; - iarg += 1; + iarg += 1; } - + } else { iarg = 4; damping_model = VISCOELASTIC; @@ -310,7 +310,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : break; } else if (strcmp(arg[iarg],"limit_damping") == 0) { limit_damping = 1; - iarg += 1; + iarg += 1; } else { error->all(FLERR, "Illegal fix wall/gran command"); } @@ -530,7 +530,7 @@ void FixWallGran::init() roll_history_index += 1; twist_history_index += 1; } - + if (damping_model == TSUJI) { double cor = normal_coeffs[1]; normal_coeffs[1] = 1.2728-4.2783*cor+11.087*pow(cor,2)-22.348*pow(cor,3)+ @@ -1020,7 +1020,7 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz, else polyhertz = sqrt((radius-r)*radius*rwall/(rwall+radius)); ccel *= polyhertz; if (limit_damping && (ccel < 0.0)) ccel = 0.0; - + // relative velocities vtr1 = vt1 - (dz*wr2-dy*wr3); @@ -1305,7 +1305,7 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz, history[thist2] -= rsht*nz; // also rescale to preserve magnitude - prjmag = sqrt(history[thist0]*history[thist0] + + prjmag = sqrt(history[thist0]*history[thist0] + history[thist1]*history[thist1] + history[thist2]*history[thist2]); if (prjmag > 0) scalefac = shrmag/prjmag; else scalefac = 0; diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp index 0f9682a133..ac26ffbf1b 100644 --- a/src/GRANULAR/pair_gran_hooke_history.cpp +++ b/src/GRANULAR/pair_gran_hooke_history.cpp @@ -370,7 +370,7 @@ void PairGranHookeHistory::settings(int narg, char **arg) xmu = utils::numeric(FLERR,arg[4],false,lmp); dampflag = utils::inumeric(FLERR,arg[5],false,lmp); if (dampflag == 0) gammat = 0.0; - + limit_damping = 0; if (narg == 7) { if (strcmp(arg[6], "limit_damping") == 0) limit_damping = 1; diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index faeaaacd88..52971f920f 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -80,7 +80,7 @@ PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp) onerad_frozen = nullptr; maxrad_dynamic = nullptr; maxrad_frozen = nullptr; - + limit_damping = nullptr; history_transfer_factors = nullptr; @@ -796,7 +796,7 @@ void PairGranular::coeff(int narg, char **arg) twist_model_one = TWIST_NONE; damping_model_one = VISCOELASTIC; int ld_flag = 0; - + int iarg = 2; while (iarg < narg) { if (strcmp(arg[iarg], "hooke") == 0) { @@ -971,7 +971,7 @@ void PairGranular::coeff(int narg, char **arg) cutoff_one = utils::numeric(FLERR,arg[iarg+1],false,lmp); iarg += 2; } else if (strcmp(arg[iarg], "limit_damping") == 0) { - ld_flag = 1; + ld_flag = 1; iarg += 1; } else error->all(FLERR, "Illegal pair_coeff command"); } @@ -1047,7 +1047,7 @@ void PairGranular::coeff(int narg, char **arg) } } - + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } @@ -1320,7 +1320,7 @@ void PairGranular::write_restart(FILE *fp) fwrite(&tangential_model[i][j],sizeof(int),1,fp); fwrite(&roll_model[i][j],sizeof(int),1,fp); fwrite(&twist_model[i][j],sizeof(int),1,fp); - fwrite(&limit_damping[i][j],sizeof(int),1,fp); + fwrite(&limit_damping[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);