From 37a046cf1eadee8d20fcfbaa6b744680bdbb0f68 Mon Sep 17 00:00:00 2001 From: "Jibril B. Coulibaly" <43829860+jibril-b-coulibaly@users.noreply.github.com> Date: Wed, 14 Aug 2019 17:39:56 -0500 Subject: [PATCH 1/6] Update pair_granular.cpp Modified PairGranular::single function to return the total normal force into argument fforce. This was done for pair styles gran/* but not for the granular pari_style, resulting in the variable fforce being uninitialized. --- src/GRANULAR/pair_granular.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index b87e64a456..2813035ebb 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -1630,7 +1630,11 @@ double PairGranular::single(int i, int j, int itype, int jtype, magtortwist = -Mtcrit * signtwist; // eq 34 } } - + + // set force and return no energy + + fforce = Fntot*rinv; + // set single_extra quantities svector[0] = fs1; From cc14103f28f65a8610780f6e1f73676214e95cca Mon Sep 17 00:00:00 2001 From: "Jibril B. Coulibaly" Date: Fri, 16 Aug 2019 15:22:15 -0500 Subject: [PATCH 2/6] Bug fixes in granular pair style: - correct formula for tangent forces in style with no history in compute() and in single() functions - remove tangent history update in the single() function - implement correct output for tangent, normal and rolling forces in single() function - correct typos in documentation --- doc/src/pair_granular.txt | 8 +++---- src/GRANULAR/pair_granular.cpp | 42 +++++++++++++++------------------- 2 files changed, 23 insertions(+), 27 deletions(-) diff --git a/doc/src/pair_granular.txt b/doc/src/pair_granular.txt index f16cd9fe0b..ccfe805b67 100644 --- a/doc/src/pair_granular.txt +++ b/doc/src/pair_granular.txt @@ -100,7 +100,7 @@ on particle {i} due to contact with particle {j} is given by: \mathbf\{F\}_\{ne, Hooke\} = k_N \delta_\{ij\} \mathbf\{n\} \end\{equation\} -Where \(\delta = R_i + R_j - \|\mathbf\{r\}_\{ij\}\|\) is the particle +Where \(\delta_\{ij\} = R_i + R_j - \|\mathbf\{r\}_\{ij\}\|\) is the particle overlap, \(R_i, R_j\) are the particle radii, \(\mathbf\{r\}_\{ij\} = \mathbf\{r\}_i - \mathbf\{r\}_j\) is the vector separating the two particle centers (note the i-j ordering so that \(F_\{ne\}\) is @@ -411,8 +411,8 @@ option by an additional factor of {a}, the radius of the contact region. The tan \mathbf\{F\}_t = -min(\mu_t F_\{n0\}, \|-k_t a \mathbf\{\xi\} + \mathbf\{F\}_\mathrm\{t,damp\}\|) \mathbf\{t\} \end\{equation\} -Here, {a} is the radius of the contact region, given by \(a = \delta -R\) for all normal contact models, except for {jkr}, where it is given +Here, {a} is the radius of the contact region, given by \(a =\sqrt\{R\delta\}\) + for all normal contact models, except for {jkr}, where it is given implicitly by \(\delta = a^2/R - 2\sqrt\{\pi \gamma a/E\}\), see discussion above. To match the Mindlin solution, one should set \(k_t = 8G\), where \(G\) is the shear modulus, related to Young's modulus @@ -680,7 +680,7 @@ The single() function of these pair styles returns 0.0 for the energy of a pairwise interaction, since energy is not conserved in these dissipative potentials. It also returns only the normal component of the pairwise interaction force. However, the single() function also -calculates 10 extra pairwise quantities. The first 3 are the +calculates 12 extra pairwise quantities. The first 3 are the components of the tangential force between particles I and J, acting on particle I. The 4th is the magnitude of this tangential force. The next 3 (5-7) are the components of the rolling torque acting on diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 2813035ebb..334c6a471e 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -391,6 +391,7 @@ void PairGranular::compute(int eflag, int vflag) } else { Fncrit = fabs(Fntot); } + Fscrit = tangential_coeffs[itype][jtype][2] * Fncrit; //------------------------------ // tangential forces @@ -446,7 +447,6 @@ void PairGranular::compute(int eflag, int vflag) fs3 = -k_tangential*history[2] - damp_tangential*vtr3; // rescale frictional displacements and forces if needed - Fscrit = tangential_coeffs[itype][jtype][2] * Fncrit; fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); if (fs > Fscrit) { shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + @@ -464,8 +464,8 @@ void PairGranular::compute(int eflag, int vflag) } 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; // From documentation: F_{t,damp} = - \eta_t v_{t,rel}, no need for extra `meff` + if (vrel != 0.0) Ft = MIN(Fscrit,fs) / vrel; // From documentation: critical force `Fscrit` used, not elastic normal force `Fne` else Ft = 0.0; fs1 = -Ft*vtr1; fs2 = -Ft*vtr2; @@ -635,7 +635,7 @@ void PairGranular::compute(int eflag, int vflag) torque[j][2] -= torroll3; } } - if (evflag) ev_tally_xyz(i,j,nlocal,0, + if (evflag) ev_tally_xyz(i,j,nlocal,0,//Should `newton_pair` passed instead of 0 ? 0.0,0.0,fx,fy,fz,delx,dely,delz); } } @@ -1451,11 +1451,13 @@ double PairGranular::single(int i, int j, int itype, int jtype, } if (damping_model[itype][jtype] == VELOCITY) { - damp_normal = normal_coeffs[itype][jtype][1]; + damp_normal = 1; + } else if (damping_model[itype][jtype] == MASS_VELOCITY) { + damp_normal = meff; } else if (damping_model[itype][jtype] == VISCOELASTIC) { - damp_normal = normal_coeffs[itype][jtype][1]*a*meff; + damp_normal = a*meff; } else if (damping_model[itype][jtype] == TSUJI) { - damp_normal = normal_coeffs[itype][jtype][1]*sqrt(meff*knfac); + damp_normal = sqrt(meff*knfac); } damp_normal_prefactor = normal_coeffs[itype][jtype][1]*damp_normal; @@ -1473,6 +1475,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, if (neighprev >= jnum) neighprev = 0; if (jlist[neighprev] == j) break; } + // the `history` pointer must not be modified here in single() function. already calculated in the compute() function. If modified here it changes the pair forces that have friction/twisting/rolling and history effects ! history = &allhistory[size_history*neighprev]; } @@ -1506,6 +1509,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, } else { Fncrit = fabs(Fntot); } + Fscrit = tangential_coeffs[itype][jtype][2] * Fncrit; //------------------------------ // tangential forces @@ -1518,13 +1522,6 @@ double PairGranular::single(int i, int j, int itype, int jtype, k_tangential *= a; } else if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) { k_tangential *= a; - // on unloading, rescale the shear displacements - if (a < history[3]) { - double factor = a/history[3]; - history[0] *= factor; - history[1] *= factor; - history[2] *= factor; - } } shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + @@ -1535,28 +1532,26 @@ double PairGranular::single(int i, int j, int itype, int jtype, fs2 = -k_tangential*history[1] - damp_tangential*vtr2; fs3 = -k_tangential*history[2] - damp_tangential*vtr3; - // rescale frictional displacements and forces if needed - Fscrit = tangential_coeffs[itype][jtype][2] * Fncrit; + // rescale frictional forces if needed fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); if (fs > Fscrit) { if (shrmag != 0.0) { - history[0] = -1.0/k_tangential*(Fscrit*fs1/fs + damp_tangential*vtr1); - history[1] = -1.0/k_tangential*(Fscrit*fs2/fs + damp_tangential*vtr2); - history[2] = -1.0/k_tangential*(Fscrit*fs3/fs + damp_tangential*vtr3); fs1 *= Fscrit/fs; fs2 *= Fscrit/fs; fs3 *= Fscrit/fs; - } else fs1 = fs2 = fs3 = 0.0; + fs *= Fscrit/fs; // saves the correct value of `fs` to svector + } else fs1 = fs2 = fs3 = fs = 0.0; // saves the correct of `fs` value to svector } // classic pair gran/hooke (no history) } else { - 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; fs3 = -Ft*vtr3; + fs = Ft*vrel; // saves the correct value of `fs` to svector } //**************************************** @@ -1601,7 +1596,8 @@ double PairGranular::single(int i, int j, int itype, int jtype, fr1 *= Frcrit/fr; fr2 *= Frcrit/fr; fr3 *= Frcrit/fr; - } else fr1 = fr2 = fr3 = 0.0; + fr *= Frcrit/fr; // saves the correct value of `fr` to svector + } else fr1 = fr2 = fr3 = fr = 0.0; // saves the correct value of `fr` to svector } } From a5acf1655bd05de4d13921d25fb647f1d5d29ab3 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 16 Aug 2019 17:30:37 -0400 Subject: [PATCH 3/6] resolve small formatting glitch Text blocks must all be flush on the left side or else sphinx gets confused since indenting is part of the syntax. --- doc/src/pair_granular.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/pair_granular.txt b/doc/src/pair_granular.txt index ccfe805b67..9fcc4dbe48 100644 --- a/doc/src/pair_granular.txt +++ b/doc/src/pair_granular.txt @@ -412,7 +412,7 @@ option by an additional factor of {a}, the radius of the contact region. The tan \end\{equation\} Here, {a} is the radius of the contact region, given by \(a =\sqrt\{R\delta\}\) - for all normal contact models, except for {jkr}, where it is given +for all normal contact models, except for {jkr}, where it is given implicitly by \(\delta = a^2/R - 2\sqrt\{\pi \gamma a/E\}\), see discussion above. To match the Mindlin solution, one should set \(k_t = 8G\), where \(G\) is the shear modulus, related to Young's modulus From bfacab7be7565f5ac8f21b21cf9f18ba3e47e3cc Mon Sep 17 00:00:00 2001 From: "Jibril B. Coulibaly" Date: Wed, 25 Sep 2019 08:55:33 -0500 Subject: [PATCH 4/6] remove non-compliant comments --- src/GRANULAR/pair_granular.cpp | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 334c6a471e..0b1a6eb09a 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -464,8 +464,8 @@ void PairGranular::compute(int eflag, int vflag) } else fs1 = fs2 = fs3 = 0.0; } } else { // classic pair gran/hooke (no history) - fs = damp_tangential*vrel; // From documentation: F_{t,damp} = - \eta_t v_{t,rel}, no need for extra `meff` - if (vrel != 0.0) Ft = MIN(Fscrit,fs) / vrel; // From documentation: critical force `Fscrit` used, not elastic normal force `Fne` + fs = damp_tangential*vrel; + if (vrel != 0.0) Ft = MIN(Fscrit,fs) / vrel; else Ft = 0.0; fs1 = -Ft*vtr1; fs2 = -Ft*vtr2; @@ -635,7 +635,7 @@ void PairGranular::compute(int eflag, int vflag) torque[j][2] -= torroll3; } } - if (evflag) ev_tally_xyz(i,j,nlocal,0,//Should `newton_pair` passed instead of 0 ? + if (evflag) ev_tally_xyz(i,j,nlocal,0, 0.0,0.0,fx,fy,fz,delx,dely,delz); } } @@ -1475,7 +1475,6 @@ double PairGranular::single(int i, int j, int itype, int jtype, if (neighprev >= jnum) neighprev = 0; if (jlist[neighprev] == j) break; } - // the `history` pointer must not be modified here in single() function. already calculated in the compute() function. If modified here it changes the pair forces that have friction/twisting/rolling and history effects ! history = &allhistory[size_history*neighprev]; } @@ -1539,8 +1538,8 @@ double PairGranular::single(int i, int j, int itype, int jtype, fs1 *= Fscrit/fs; fs2 *= Fscrit/fs; fs3 *= Fscrit/fs; - fs *= Fscrit/fs; // saves the correct value of `fs` to svector - } else fs1 = fs2 = fs3 = fs = 0.0; // saves the correct of `fs` value to svector + fs *= Fscrit/fs; + } else fs1 = fs2 = fs3 = fs = 0.0; } // classic pair gran/hooke (no history) @@ -1551,7 +1550,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, fs1 = -Ft*vtr1; fs2 = -Ft*vtr2; fs3 = -Ft*vtr3; - fs = Ft*vrel; // saves the correct value of `fs` to svector + fs = Ft*vrel; } //**************************************** @@ -1596,8 +1595,8 @@ double PairGranular::single(int i, int j, int itype, int jtype, fr1 *= Frcrit/fr; fr2 *= Frcrit/fr; fr3 *= Frcrit/fr; - fr *= Frcrit/fr; // saves the correct value of `fr` to svector - } else fr1 = fr2 = fr3 = fr = 0.0; // saves the correct value of `fr` to svector + fr *= Frcrit/fr; + } else fr1 = fr2 = fr3 = fr = 0.0; } } From 12803b7dcfd8f04cd604e0ce1f6696391f34fa5e Mon Sep 17 00:00:00 2001 From: "Jibril B. Coulibaly" Date: Thu, 26 Sep 2019 16:01:24 -0500 Subject: [PATCH 5/6] add newton_pair flag to ev_tally_xyz function --- src/GRANULAR/pair_granular.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 0b1a6eb09a..43bcec5181 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -635,7 +635,7 @@ void PairGranular::compute(int eflag, int vflag) torque[j][2] -= torroll3; } } - if (evflag) ev_tally_xyz(i,j,nlocal,0, + if (evflag) ev_tally_xyz(i,j,nlocal,force->newton_pair, 0.0,0.0,fx,fy,fz,delx,dely,delz); } } From 747bfa31f26c5c10e549f36d46451989a9c5ccae Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 3 Oct 2019 18:31:42 +0200 Subject: [PATCH 6/6] apply changes suggested by Andrew Santos --- doc/src/pair_granular.txt | 14 +++++++------- src/GRANULAR/pair_granular.cpp | 27 ++++++++------------------- 2 files changed, 15 insertions(+), 26 deletions(-) diff --git a/doc/src/pair_granular.txt b/doc/src/pair_granular.txt index 9fcc4dbe48..718315b4fe 100644 --- a/doc/src/pair_granular.txt +++ b/doc/src/pair_granular.txt @@ -177,7 +177,7 @@ following general form: \end\{equation\} Here, \(\mathbf\{v\}_\{n,rel\} = (\mathbf\{v\}_j - \mathbf\{v\}_i) -\cdot \mathbf\{n\}\) is the component of relative velocity along +\cdot \mathbf\{n\} \mathbf\{n\}\) is the component of relative velocity along \(\mathbf\{n\}\). The optional {damping} keyword to the {pair_coeff} command followed by @@ -299,8 +299,8 @@ the normal damping \(\eta_n\) (see above): \eta_t = -x_\{\gamma,t\} \eta_n \end\{equation\} -The normal damping prefactor \(\eta_n\) is determined by the choice of -the {damping} keyword, as discussed above. Thus, the {damping} +The normal damping prefactor \(\eta_n\) is determined by the choice +of the {damping} keyword, as discussed above. Thus, the {damping} keyword also affects the tangential damping. The parameter \(x_\{\gamma,t\}\) is a scaling coefficient. Several works in the literature use \(x_\{\gamma,t\} = 1\) ("Marshall"_#Marshall2009, @@ -308,10 +308,10 @@ literature use \(x_\{\gamma,t\} = 1\) ("Marshall"_#Marshall2009, tangential velocity at the point of contact is given by \(\mathbf\{v\}_\{t, rel\} = \mathbf\{v\}_\{t\} - (R_i\Omega_i + R_j\Omega_j) \times \mathbf\{n\}\), where \(\mathbf\{v\}_\{t\} = -\mathbf\{v\}_r - \mathbf\{v\}_r\cdot\mathbf\{n\}\), \(\mathbf\{v\}_r = -\mathbf\{v\}_j - \mathbf\{v\}_i\). The direction of the applied force -is \(\mathbf\{t\} = -\mathbf\{v_\{t,rel\}\}/\|\mathbf\{v_\{t,rel\}\}\|\). +\mathbf\{v\}_r - \mathbf\{v\}_r\cdot\mathbf\{n\}\{n\}\), +\(\mathbf\{v\}_r = \mathbf\{v\}_j - \mathbf\{v\}_i\). +The direction of the applied force is \(\mathbf\{t\} = +\mathbf\{v_\{t,rel\}\}/\|\mathbf\{v_\{t,rel\}\}\|\) . The normal force value \(F_\{n0\}\) used to compute the critical force depends on the form of the contact model. For non-cohesive models diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 43bcec5181..54f6b77d2f 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -305,7 +305,8 @@ void PairGranular::compute(int eflag, int vflag) delta = radsum - r; dR = delta*Reff; - if (normal_model[itype][jtype] == JKR) { + + if (normal_model[itype][jtype] == JKR) { touch[jj] = 1; R2=Reff*Reff; coh = normal_coeffs[itype][jtype][3]; @@ -1374,22 +1375,6 @@ double PairGranular::single(int i, int j, int itype, int jtype, vn2 = ny*vnnr; vn3 = nz*vnnr; - double *rmass = atom->rmass; - int *mask = atom->mask; - mi = rmass[i]; - mj = rmass[j]; - if (fix_rigid) { - if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; - if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; - } - - meff = mi*mj / (mi+mj); - if (mask[i] & freeze_group_bit) meff = mj; - if (mask[j] & freeze_group_bit) meff = mi; - - delta = radsum - r; - dR = delta*Reff; - // tangential component vt1 = vr1 - vn1; @@ -1407,6 +1392,9 @@ double PairGranular::single(int i, int j, int itype, int jtype, // 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) { @@ -1557,7 +1545,8 @@ double PairGranular::single(int i, int j, int itype, int jtype, // rolling resistance //**************************************** - if (roll_model[itype][jtype] != ROLL_NONE) { + 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]; @@ -1623,7 +1612,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, Mtcrit = mu_twist*Fncrit; // critical torque (eq 44) if (fabs(magtortwist) > Mtcrit) { magtortwist = -Mtcrit * signtwist; // eq 34 - } + } else magtortwist = 0.0; } // set force and return no energy