From ef803be08e65999b2b1c657d32ec6187ec774569 Mon Sep 17 00:00:00 2001 From: dsbolin Date: Tue, 15 Jan 2019 10:18:46 -0700 Subject: [PATCH] Moved all model option syntax for pair granular to pair coeff command; added global cutoff option for pair style granular command; initial write-up of documentation. --- doc/src/pair_granular.txt | 332 +++++++++++++++++++++++ src/GRANULAR/pair_granular.cpp | 392 ++++++++------------------- src/GRANULAR/pair_granular.h | 13 +- src/GRANULAR/pair_granular_multi.cpp | 379 ++++++-------------------- src/GRANULAR/pair_granular_multi.h | 17 +- 5 files changed, 538 insertions(+), 595 deletions(-) create mode 100644 doc/src/pair_granular.txt diff --git a/doc/src/pair_granular.txt b/doc/src/pair_granular.txt new file mode 100644 index 0000000000..19c16cc6ea --- /dev/null +++ b/doc/src/pair_granular.txt @@ -0,0 +1,332 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Commands_all.html) + +:line + +pair_style granular command :h3 +pair_style granular/multi command :h3 + +[Syntax:] + +pair_style style cutoff :pre + +style = {granular} or {granular/multi} :ulb,l +cutoff = global cutoff (optional). See discussion below. + +[Examples:] + +pair_style granular +pair coeff 1 1 hertz 1000.0 50.0 tangential mindlin 800.0 50.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall +pair coeff 2 2 hertz 200.0 20.0 tangential mindlin 300.0 50.0 0.1 rolling sds 200.0 100.0 0.1 twisting marshall + +pair_style granular/multi +pair coeff 1 1 hertz 1000.0 50.0 tangential mindlin 800.0 50.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall +pair coeff 2 2 dmt 1000.0 50.0 800.0 10.0 tangential mindlin 800.0 50.0 0.1 roll sds 500.0 200.0 0.1 twisting marshall +pair coeff 1 2 dmt 1000.0 50.0 800.0 10.0 tangential mindlin 800.0 50.0 0.1 roll sds 500.0 200.0 0.1 twisting marshall + + +[Description:] + +The {granular} styles support a variety of options for the normal, tangential, rolling and twisting +forces resulting from contact between two granular particles. The computed force depends on the combination +of choices for these models. + +All model options and parameters are entered in the "pair_coeff"_pair_coeff.html command, as described below. +Unlike e.g. "pair gran/hooke"_pair_gran.html, coefficient values are not global, but can be set to different values for +various combinations of particle types, as determined by the "pair_coeff"_pair_coeff.html command. +In the case of {granular}, coefficients +can vary between particle types, but model choices cannot. For instance, in the first +example above, the stiffness, damping, and tangential friction are different for type 1 - type 1 and type 2 - type 2 interactions, but +both 1-1 and 2-2 interactions must have the same model form, hence all keywords are identical between the two types. Cross-coefficients +for 1-2 interactions for the case of the {hertz} model above are set via simple geometric mixing rules. The {granular/multi} +style removes this restriction at a small cost in computational efficiency, so that different particle types +can potentially interact via different model forms. As shown in the second example, +1-1 interactions are based on a Hertzian contact model and 2-2 interactions are based on a {dmt} model (see below). +In the case that 1-1 and 2-2 interactions have different model forms, mixing of coefficients cannot be +determined, so 1-2 interactions must be explicitly defined via the pair coeff command, otherwise an error results. + +The first required keyword for the pair coeff command is the normal contact model. Currently supported options and +the required arguments are: + +{hooke} : k_n, damping +{hertz} : k_n, damping +{hertz/material} : E, damping, G +{dmt} : E, damping, G, cohesion +{jkr} : E, damping, G, cohesion + +Here, k_n is spring stiffness, damping is a damping constant or a coefficient of restitution, depending on +the choice of damping model (see below), E and G are Young's modulus and shear modulus, in units of pressure, +and cohesion is a surface energy density, in units of energy/length^2. + +For the {hooke} model, the normal component of force is given by: +:c,image(Eqs/hooke_normal.jpg) + +For {hertz}, the normal force is given by: +:c,image{Eqs/hertz_normal.jpg} + +For both [hooke] and [hertz], stiffness for unspecified cross-terms is given by simple geometric mixing +(e.g. if stiffness is specified for type 1 and type 2 particles as k_n_1 and k_n_2, respectively, +type 1 - type 2 contacts use a stiffness given by k_n_{12} = sqrt(k_n_1*k_n_2)) + +For {hertz/material}, the form is the same as above, but coefficients are computed differently, and mixing follows +a different rule based on shear modulus: +:c,image{Eqs/hertz_material_normal.jpg} + +For {dmt}, the normal force is given by: +:c,image{Eqs/dmt_normal.jpg} + +Where gamma is cohesion. + +For {jkr}, the normal force is given by: +:c,image{Eqs/jkr_normal.jpg} + +The same mixing rule for stiffness as for {hertz/material} is used by both the {dmt} and {jkr} models. + +The tangential contact model must also be specified, which follows +the required {tangential} keyword. Currently supported options +and their required arguments are: + +{no_history}: k_t, tangential_damping, friction coefficient +{mindlin}: k_t, tangential_damping, friction coefficient + +For {no_history}, the tangential force is computed according to + +The total force on a particle is the sum of the normal and tangential forces from all interactions. The tangential +force also induces a torque on both particles in a contacting pair. Additionally, rolling and twisting friction +models can also be applied, which may induce additional torques (but no force). The following options are +supported for the rolling friction model + + +The first required keyword +in the pair coeff command is the choice +of normal force contact model, for which current opitons are {hooke}, {hertz} + +The {gran} styles use the following formulas for the frictional force +between two granular particles, as described in +"(Brilliantov)"_#Brilliantov, "(Silbert)"_#Silbert, and +"(Zhang)"_#Zhang3, when the distance r between two particles of radii +Ri and Rj is less than their contact distance d = Ri + Rj. There is +no force between the particles when r > d. + +The two Hookean styles use this formula: + +:c,image(Eqs/pair_gran_hooke.jpg) + +The Hertzian style uses this formula: + +:c,image(Eqs/pair_gran_hertz.jpg) + +In both equations the first parenthesized term is the normal force +between the two particles and the second parenthesized term is the +tangential force. The normal force has 2 terms, a contact force and a +damping force. The tangential force also has 2 terms: a shear force +and a damping force. The shear force is a "history" effect that +accounts for the tangential displacement between the particles for the +duration of the time they are in contact. This term is included in +pair styles {hooke/history} and {hertz/history}, but is not included +in pair style {hooke}. The tangential damping force term is included +in all three pair styles if {dampflag} is set to 1; it is not included +if {dampflag} is set to 0. + +The other quantities in the equations are as follows: + +delta = d - r = overlap distance of 2 particles +Kn = elastic constant for normal contact +Kt = elastic constant for tangential contact +gamma_n = viscoelastic damping constant for normal contact +gamma_t = viscoelastic damping constant for tangential contact +m_eff = Mi Mj / (Mi + Mj) = effective mass of 2 particles of mass Mi and Mj +Delta St = tangential displacement vector between 2 particles \ + which is truncated to satisfy a frictional yield criterion +n_ij = unit vector along the line connecting the centers of the 2 particles +Vn = normal component of the relative velocity of the 2 particles +Vt = tangential component of the relative velocity of the 2 particles :ul + +The Kn, Kt, gamma_n, and gamma_t coefficients are specified as +parameters to the pair_style command. If a NULL is used for Kt, then +a default value is used where Kt = 2/7 Kn. If a NULL is used for +gamma_t, then a default value is used where gamma_t = 1/2 gamma_n. + +The interpretation and units for these 4 coefficients are different in +the Hookean versus Hertzian equations. + +The Hookean model is one where the normal push-back force for two +overlapping particles is a linear function of the overlap distance. +Thus the specified Kn is in units of (force/distance). Note that this +push-back force is independent of absolute particle size (in the +monodisperse case) and of the relative sizes of the two particles (in +the polydisperse case). This model also applies to the other terms in +the force equation so that the specified gamma_n is in units of +(1/time), Kt is in units of (force/distance), and gamma_t is in units +of (1/time). + +The Hertzian model is one where the normal push-back force for two +overlapping particles is proportional to the area of overlap of the +two particles, and is thus a non-linear function of overlap distance. +Thus Kn has units of force per area and is thus specified in units of +(pressure). The effects of absolute particle size (monodispersity) +and relative size (polydispersity) are captured in the radii-dependent +pre-factors. When these pre-factors are carried through to the other +terms in the force equation it means that the specified gamma_n is in +units of (1/(time*distance)), Kt is in units of (pressure), and +gamma_t is in units of (1/(time*distance)). + +Note that in the Hookean case, Kn can be thought of as a linear spring +constant with units of force/distance. In the Hertzian case, Kn is +like a non-linear spring constant with units of force/area or +pressure, and as shown in the "(Zhang)"_#Zhang3 paper, Kn = 4G / +(3(1-nu)) where nu = the Poisson ratio, G = shear modulus = E / +(2(1+nu)), and E = Young's modulus. Similarly, Kt = 4G / (2-nu). +(NOTE: in an earlier version of the manual, we incorrectly stated that +Kt = 8G / (2-nu).) + +Thus in the Hertzian case Kn and Kt can be set to values that +corresponds to properties of the material being modeled. This is also +true in the Hookean case, except that a spring constant must be chosen +that is appropriate for the absolute size of particles in the model. +Since relative particle sizes are not accounted for, the Hookean +styles may not be a suitable model for polydisperse systems. + +NOTE: In versions of LAMMPS before 9Jan09, the equation for Hertzian +interactions did not include the sqrt(RiRj/Ri+Rj) term and thus was +not as accurate for polydisperse systems. For monodisperse systems, +sqrt(RiRj/Ri+Rj) is a constant factor that effectively scales all 4 +coefficients: Kn, Kt, gamma_n, gamma_t. Thus you can set the values +of these 4 coefficients appropriately in the current code to reproduce +the results of a previous Hertzian monodisperse calculation. For +example, for the common case of a monodisperse system with particles +of diameter 1, all 4 of these coefficients should now be set 2x larger +than they were previously. + +Xmu is also specified in the pair_style command and is the upper limit +of the tangential force through the Coulomb criterion Ft = xmu*Fn, +where Ft and Fn are the total tangential and normal force components +in the formulas above. Thus in the Hookean case, the tangential force +between 2 particles grows according to a tangential spring and +dash-pot model until Ft/Fn = xmu and is then held at Ft = Fn*xmu until +the particles lose contact. In the Hertzian case, a similar analogy +holds, though the spring is no longer linear. + +NOTE: Normally, xmu should be specified as a fractional value between +0.0 and 1.0, however LAMMPS allows large values (up to 1.0e4) to allow +for modeling of systems which can sustain very large tangential +forces. + +The effective mass {m_eff} is given by the formula above for two +isolated particles. If either particle is part of a rigid body, its +mass is replaced by the mass of the rigid body in the formula above. +This is determined by searching for a "fix rigid"_fix_rigid.html +command (or its variants). + +For granular styles there are no additional coefficients to set for +each pair of atom types via the "pair_coeff"_pair_coeff.html command. +All settings are global and are made via the pair_style command. +However you must still use the "pair_coeff"_pair_coeff.html for all +pairs of granular atom types. For example the command + +pair_coeff * * :pre + +should be used if all atoms in the simulation interact via a granular +potential (i.e. one of the pair styles above is used). If a granular +potential is used as a sub-style of "pair_style +hybrid"_pair_hybrid.html, then specific atom types can be used in the +pair_coeff command to determine which atoms interact via a granular +potential. + +:line + +Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are +functionally the same as the corresponding style without the suffix. +They have been optimized to run faster, depending on your available +hardware, as discussed on the "Speed packages"_Speed_packages.html doc +page. The accelerated styles take the same arguments and should +produce the same results, except for round-off and precision issues. + +These accelerated styles are part of the GPU, USER-INTEL, KOKKOS, +USER-OMP and OPT packages, respectively. They are only enabled if +LAMMPS was built with those packages. See the "Build +package"_Build_package.html doc page for more info. + +You can specify the accelerated styles explicitly in your input script +by including their suffix, or you can use the "-suffix command-line +switch"_Run_options.html when you invoke LAMMPS, or you can use the +"suffix"_suffix.html command in your input script. + +See the "Speed packages"_Speed_packages.html doc page for more +instructions on how to use the accelerated styles effectively. + +:line + +[Mixing, shift, table, tail correction, restart, rRESPA info]: + +The "pair_modify"_pair_modify.html mix, shift, table, and tail options +are not relevant for granular pair styles. + +These pair styles write their information to "binary restart +files"_restart.html, so a pair_style command does not need to be +specified in an input script that reads a restart file. + +These pair styles can only be used via the {pair} keyword of the +"run_style respa"_run_style.html command. They do not support the +{inner}, {middle}, {outer} keywords. + +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 +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 relative velocity in the +normal direction (along the line joining the 2 sphere centers). The +last 3 (8-10) the components of the relative velocity in the +tangential direction. + +These extra quantities can be accessed by the "compute +pair/local"_compute_pair_local.html command, as {p1}, {p2}, ..., +{p10}. + +:line + +[Restrictions:] + +All the granular pair styles are part of the GRANULAR package. It is +only enabled if LAMMPS was built with that package. See the "Build +package"_Build_package.html doc page for more info. + +These pair styles require that atoms store torque and angular velocity +(omega) as defined by the "atom_style"_atom_style.html. They also +require a per-particle radius is stored. The {sphere} atom style does +all of this. + +This pair style requires you to use the "comm_modify vel +yes"_comm_modify.html command so that velocities are stored by ghost +atoms. + +These pair styles will not restart exactly when using the +"read_restart"_read_restart.html command, though they should provide +statistically similar results. This is because the forces they +compute depend on atom velocities. See the +"read_restart"_read_restart.html command for more details. + +[Related commands:] + +"pair_coeff"_pair_coeff.html + +[Default:] none + +:line + +:link(Brilliantov) +[(Brilliantov)] Brilliantov, Spahn, Hertzsch, Poschel, Phys Rev E, 53, +p 5382-5392 (1996). + +:link(Silbert) +[(Silbert)] Silbert, Ertas, Grest, Halsey, Levine, Plimpton, Phys Rev +E, 64, p 051302 (2001). + +:link(Zhang3) +[(Zhang)] Zhang and Makse, Phys Rev E, 72, p 011301 (2005). diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index c75c80fea5..9fe4bd8415 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -1360,205 +1360,15 @@ void PairGranular::allocate() void PairGranular::settings(int narg, char **arg) { - if (narg < 5) error->all(FLERR,"Illegal pair_style command"); - - int iarg = 0; - - //Some defaults - normal = HERTZ; - damping = VISCOELASTIC; - tangential = TANGENTIAL_MINDLIN; - roll = ROLL_NONE; - twist = TWIST_NONE; - - tangential_history = 1; + if (narg == 1){ + cutoff_global = force->numeric(FLERR,arg[0]) + } + else{ + cutoff_global = -1; //Will be set based on particle sizes, model choice + } + tangential_history = 0; roll_history = twist_history = 0; - - int normal_set, tangential_set; - normal_set = tangential_set = 0; - - while (iarg < narg){ - if (strcmp(arg[iarg], "hooke") == 0){ - if (iarg + 2 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for Hooke option"); - normal = HOOKE; - memory->create(normal_coeffs_global, 2, "pair:normal_coeffs_global"); - normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1]); //kn - normal_coeffs_global[1] = force->numeric(FLERR,arg[iarg+2]); //damping - iarg += 3; - } - else if (strcmp(arg[iarg], "hertz") == 0){ - int num_coeffs = 2; - if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for Hertz option"); - normal = HERTZ; - normal_set = 1; - memory->create(normal_coeffs_global, num_coeffs, "pair:normal_coeffs_global"); - normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1]); //kn - normal_coeffs_global[1] = force->numeric(FLERR,arg[iarg+2]); //damping - iarg += num_coeffs+1; - } - else if (strcmp(arg[iarg], "hertz/material") == 0){ - int num_coeffs = 3; - if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for Hertz option"); - normal = HERTZ_MATERIAL; - normal_set = 1; - memory->create(normal_coeffs_global, num_coeffs, "pair:normal_coeffs_global"); - normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //E (Young's modulus) - normal_coeffs_global[1] = force->numeric(FLERR,arg[iarg+2]); //damping - normal_coeffs_global[2] = force->numeric(FLERR,arg[iarg+3]); //G (shear modulus) - iarg += num_coeffs+1; - } - else if (strcmp(arg[iarg], "dmt") == 0){ - if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for Hertz option"); - normal = DMT; - normal_set = 1; - memory->create(normal_coeffs_global, 4, "pair:normal_coeffs_global"); - normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //4/3 E - normal_coeffs_global[1] = force->numeric(FLERR,arg[iarg+2]); //damping - normal_coeffs_global[2] = force->numeric(FLERR,arg[iarg+3]); //G - normal_coeffs_global[3] = force->numeric(FLERR,arg[iarg+3]); //cohesion - iarg += 5; - } - else if (strcmp(arg[iarg], "jkr") == 0){ - if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for JKR option"); - beyond_contact = 1; - normal = JKR; - normal_set = 1; - memory->create(normal_coeffs_global, 4, "pair:normal_coeffs_global"); - normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1]); //E - normal_coeffs_global[1] = force->numeric(FLERR,arg[iarg+2]); //damping - normal_coeffs_global[2] = force->numeric(FLERR,arg[iarg+3]); //G - normal_coeffs_global[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion - iarg += 5; - } - else if (strcmp(arg[iarg], "damp_velocity") == 0){ - damping = VELOCITY; - iarg += 1; - } - else if (strcmp(arg[iarg], "damp_viscoelastic") == 0){ - damping = VISCOELASTIC; - iarg += 1; - } - else if (strcmp(arg[iarg], "damp_tsuji") == 0){ - damping = TSUJI; - iarg += 1; - } - else if (strstr(arg[iarg], "tangential") != NULL){ - if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for tangential model"); - if (strstr(arg[iarg+1], "nohistory") != NULL){ - tangential = TANGENTIAL_NOHISTORY; - tangential_set = 1; - } - else if (strstr(arg[iarg+1], "mindlin") != NULL){ - tangential = TANGENTIAL_MINDLIN; - tangential_history = 1; - tangential_set = 1; - } - else{ - error->all(FLERR, "Illegal pair_style command, unrecognized sliding friction model"); - } - memory->create(tangential_coeffs_global, 3, "pair:tangential_coeffs_global"); - tangential_coeffs_global[0] = force->numeric(FLERR,arg[iarg+2]); //kt - tangential_coeffs_global[1] = force->numeric(FLERR,arg[iarg+3]); //gammat - tangential_coeffs_global[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. - iarg += 5; - } - else if (strstr(arg[iarg], "roll") != NULL){ - if (strstr(arg[iarg+1], "none") != NULL){ - roll = ROLL_NONE; - iarg += 2; - } - else{ - if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for rolling model"); - if (strstr(arg[iarg+1], "nohistory") != NULL){ - roll = ROLL_NOHISTORY; - } - else if (strstr(arg[iarg+1], "sds") != NULL){ - roll = ROLL_SDS; - roll_history = 1; - } - else{ - error->all(FLERR, "Illegal pair_style command, unrecognized rolling friction model"); - } - memory->create(roll_coeffs_global, 3, "pair:roll_coeffs_global"); - roll_coeffs_global[0] = force->numeric(FLERR,arg[iarg+2]); //kR - roll_coeffs_global[1] = force->numeric(FLERR,arg[iarg+3]); //gammaR - roll_coeffs_global[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. - iarg += 5; - } - } - else if (strstr(arg[iarg], "twist") != NULL){ - if (strstr(arg[iarg+1], "none") != NULL){ - twist = TWIST_NONE; - iarg += 2; - } - else if (strstr(arg[iarg+1], "marshall") != NULL){ - twist = TWIST_MARSHALL; - twist_history = 1; - iarg += 2; - } - else{ - if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for twist model"); - memory->create(twist_coeffs_global, 3, "pair:twist_coeffs_global"); //To be filled later - if (strstr(arg[iarg+1], "nohistory") != NULL){ - twist = TWIST_NOHISTORY; - } - else if (strstr(arg[iarg+1], "sds") != NULL){ - twist = TWIST_SDS; - twist_history = 1; - } - else{ - error->all(FLERR, "Illegal pair_style command, unrecognized twisting friction model"); - } - memory->create(twist_coeffs_global, 3, "pair:twist_coeffs_global"); - twist_coeffs_global[0] = force->numeric(FLERR,arg[iarg+2]); //ktwist - twist_coeffs_global[1] = force->numeric(FLERR,arg[iarg+3]); //gammatwist - twist_coeffs_global[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. - iarg += 5; - } - } - else error->all(FLERR, "Illegal pair_style granular command"); - } - - //Set all i-i entries, which may be replaced by pair coeff commands - //It may also make sense to consider removing all of the above, and only - // having the option for pair_coeff to set the parameters, similar to most LAMMPS pair styles - // The reason for the current setup is to remain true to existing pair gran/hooke etc. syntax, - // where coeffs are set in the pair_style command, and a pair_coeff * * command is issued. - - allocate(); - double damp; - if (damping == TSUJI){ - double cor = normal_coeffs_global[1]; - damp = 1.2728-4.2783*cor+11.087*pow(cor,2)-22.348*pow(cor,3)+ - 27.467*pow(cor,4)-18.022*pow(cor,5)+ - 4.8218*pow(cor,6); - } - else damp = normal_coeffs_global[1]; - - for (int i = 1; i <= atom->ntypes; i++){ - if (normal_set){ - normal_coeffs[i][i][0] = normal_coeffs_global[0]; - normal_coeffs[i][i][1] = damp; - if (normal != HOOKE && normal != HERTZ){ - normal_coeffs[i][i][2] = normal_coeffs_global[2]; - } - if ((normal == JKR) || (normal == DMT)) - normal_coeffs[i][i][3] = normal_coeffs_global[3]; - } - if(tangential_set){ - for (int k = 0; k < 3; k++) - tangential_coeffs[i][i][k] = tangential_coeffs_global[k]; - } - if (roll != ROLL_NONE) - for (int k = 0; k < 3; k++) - roll_coeffs[i][i][k] = roll_coeffs_global[k]; - - if (twist != TWIST_NONE && twist != TWIST_MARSHALL) - for (int k = 0; k < 3; k++) - twist_coeffs[i][i][k] = twist_coeffs_global[k]; - - if (normal_set && tangential_set) setflag[i][i] = 1; - } + normal_set = tangential_set = damping_set = roll_set = twist_set = 0; } /* ---------------------------------------------------------------------- @@ -1567,12 +1377,6 @@ void PairGranular::settings(int narg, char **arg) void PairGranular::coeff(int narg, char **arg) { - int normal_set, damping_set, tangential_set, roll_set, twist_set; - double *normal_coeffs_local; - double *tangential_coeffs_local; - double *roll_coeffs_local; - double *twist_coeffs_local; - normal_coeffs_local = new double[4]; tangential_coeffs_local = new double[4]; roll_coeffs_local = new double[4]; @@ -1587,13 +1391,12 @@ void PairGranular::coeff(int narg, char **arg) force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi); force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi); - normal_set = damping_set = tangential_set = roll_set = twist_set = 0; - 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"); - if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be consistent"); + if (!normal_set) normal = HOOKE; + else if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be the same for all types"); normal_coeffs_local[0] = force->numeric(FLERR,arg[iarg+1]); //kn normal_coeffs_local[1] = force->numeric(FLERR,arg[iarg+2]); //damping normal_set = 1; @@ -1602,7 +1405,8 @@ void PairGranular::coeff(int narg, char **arg) else if (strcmp(arg[iarg], "hertz") == 0){ int num_coeffs = 2; if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hertz option"); - if (normal != HERTZ) if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be consistent"); + if (!normal_set) normal = HERTZ; + else if (normal_set && normal != HERTZ) if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be the same for all types"); normal_coeffs_local[0] = force->numeric(FLERR,arg[iarg+1]); //kn normal_coeffs_local[1] = force->numeric(FLERR,arg[iarg+2]); //damping normal_set = 1; @@ -1611,7 +1415,8 @@ void PairGranular::coeff(int narg, char **arg) else if (strcmp(arg[iarg], "hertz/material") == 0){ int num_coeffs = 3; if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hertz option"); - if (normal != HERTZ) if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be consistent"); + if (!normal_set) normal = HERTZ/MATERIAL; + else if (normal != HERTZ) if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be the same for all types"); normal_coeffs_local[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //E normal_coeffs_local[1] = force->numeric(FLERR,arg[iarg+2]); //damping normal_coeffs_local[2] = force->numeric(FLERR,arg[iarg+3]); //G @@ -1620,7 +1425,8 @@ void PairGranular::coeff(int narg, char **arg) } 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"); - if (normal != DMT) if (normal != HERTZ) if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be consistent"); + if (!normal_set) normal = DMT; + else if (normal != DMT) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be the same for all types"); normal_coeffs_local[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //E normal_coeffs_local[1] = force->numeric(FLERR,arg[iarg+2]); //damping normal_coeffs_local[2] = force->numeric(FLERR,arg[iarg+3]); //G @@ -1631,7 +1437,8 @@ void PairGranular::coeff(int narg, char **arg) 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; - if (normal != JKR) if (normal != HERTZ) if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be consistent"); + if (!normal_set) normal = JKR; + else if (normal != JKR) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be the same for all types"); normal_coeffs_local[0] = force->numeric(FLERR,arg[iarg+1]); //E normal_coeffs_local[1] = force->numeric(FLERR,arg[iarg+2]); //damping normal_coeffs_local[2] = force->numeric(FLERR,arg[iarg+3]); //G @@ -1639,25 +1446,32 @@ void PairGranular::coeff(int narg, char **arg) normal_set = 1; iarg += 5; } - else if (strcmp(arg[iarg], "damp_velocity") == 0){ - if (damping != VELOCITY) error->all(FLERR, "Illegal pair_coeff command, choice of damping contact model must be consistent"); - iarg += 1; - } - else if (strcmp(arg[iarg], "damp_viscoelastic") == 0){ - if (damping != VISCOELASTIC) error->all(FLERR, "Illegal pair_coeff command, choice of damping contact model must be consistent"); - iarg += 1; - } - else if (strcmp(arg[iarg], "damp_tsuji") == 0){ - if (damping != TSUJI) error->all(FLERR, "Illegal pair_coeff command, choice of damping contact model must be consistent"); - iarg += 1; - } - else if (strstr(arg[iarg], "tangential") != NULL){ - if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model"); - if (strstr(arg[iarg+1], "nohistory") != NULL){ - if (tangential != TANGENTIAL_NOHISTORY) error->all(FLERR, "Illegal pair_coeff command, choice of tangential contact model must be consistent"); + else if (strcmp(arg[iarg], "damp") == 0){ + if (iarg+1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters provided for damping model"); + if (strcmp(arg[iarg+1]), "velocity") == 0){ + if (!damping_set) damping = VELOCITY; + else if (damping != VELOCITY) error->all(FLERR, "Illegal pair_coeff command, choice of damping contact model must be the same for all types"); } - else if (strstr(arg[iarg+1], "mindlin") != NULL){ - if (tangential != TANGENTIAL_MINDLIN) error->all(FLERR, "Illegal pair_coeff command, choice of tangential contact model must be consistent");; + else if (strcmp(arg[iarg+1], "viscoelastic") == 0){ + if (!damping_set) damping = VISCOELASTIC; + else if (damping != VISCOELASTIC) error->all(FLERR, "Illegal pair_coeff command, choice of damping contact model must be the same for all types"); + } + else if (strcmp(arg[iarg+1], "tsuji") == 0){ + if (!damping_set) damping = TSUJI; + if (damping != TSUJI) error->all(FLERR, "Illegal pair_coeff command, choice of damping contact model must be the same for all types"); + } + damping_set = 1; + iarg += 1; + } + else if (strcmp(arg[iarg], "tangential") == 0){ + if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model"); + if (strcmp(arg[iarg+1], "nohistory") == 0){ + if (!tangential_set) tangential = TANGENTIAL_NOHISTORY; + else if (tangential != TANGENTIAL_NOHISTORY) error->all(FLERR, "Illegal pair_coeff command, choice of tangential contact model must be the same for all types"); + } + else if (strcmp(arg[iarg+1], "mindlin") == 0){ + if (!tangential_set) tangential = TANGENTIAL_MINDLIN; + else if (tangential != TANGENTIAL_MINDLIN) error->all(FLERR, "Illegal pair_coeff command, choice of tangential contact model must be the same for all types");; tangential_history = 1; } else{ @@ -1669,19 +1483,22 @@ void PairGranular::coeff(int narg, char **arg) tangential_coeffs_local[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. iarg += 5; } - else if (strstr(arg[iarg], "rolling") != NULL){ + else if (strcmp(arg[iarg], "rolling") == 0){ if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); - if (strstr(arg[iarg+1], "none") != NULL){ - if (roll != ROLL_NONE) error->all(FLERR, "Illegal pair_coeff command, choice of rolling friction model must be consistent"); + if (strcmp(arg[iarg+1], "none") == 0){ + if (!roll_set) roll = ROLL_NONE; + else if (roll != ROLL_NONE) error->all(FLERR, "Illegal pair_coeff command, choice of rolling friction model must be the same for all types"); iarg += 2; } else{ if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for rolling model"); - if (strstr(arg[iarg+1], "nohistory") != NULL){ - if (roll != ROLL_NOHISTORY) error->all(FLERR, "Illegal pair_coeff command, choice of rolling friction model must be consistent"); + if (strcmp(arg[iarg+1], "nohistory") == 0){ + if (!roll_set) roll = ROLL_NOHISTORY; + else if (roll != ROLL_NOHISTORY) error->all(FLERR, "Illegal pair_coeff command, choice of rolling friction model must be the same for all types"); } - else if (strstr(arg[iarg+1], "sds") != NULL){ - if (roll != ROLL_SDS) error->all(FLERR, "Illegal pair_coeff command, choice of rolling friction model must be consistent"); + else if (strcmp(arg[iarg+1], "sds") == 0){ + if (!roll_set) roll = ROLL_SDS; + else if (roll != ROLL_SDS) error->all(FLERR, "Illegal pair_coeff command, choice of rolling friction model must be the same for all types"); roll_history = 1; } else{ @@ -1694,24 +1511,28 @@ void PairGranular::coeff(int narg, char **arg) iarg += 5; } } - else if (strstr(arg[iarg], "twist") != NULL){ + else if (strcmp(arg[iarg], "twisting") == 0){ if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); - if (strstr(arg[iarg+1], "none") != NULL){ - if (twist != TWIST_NONE) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be consistent"); + if (strcmp(arg[iarg+1], "none") == 0){ + if (!twist_set) twist = TWIST_NOHISTORY; + else if (twist != TWIST_NONE) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be the same for all types"); iarg += 2; } - else if (strstr(arg[iarg+1], "marshall") != NULL){ - if (twist != TWIST_MARSHALL) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be consistent"); + else if (strcmp(arg[iarg+1], "marshall") == 0){ + if (!twist_set) twist = TWIST_MARSHALL; + else if (twist != TWIST_MARSHALL) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be the same for all types"); twist_history = 1; iarg += 2; } else{ if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for twist model"); - if (strstr(arg[iarg+1], "nohistory") != NULL){ - if (twist != TWIST_NOHISTORY) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be consistent"); + if (!twist_set) twist = TWIST_NOHISTORY; + else if (strcmp(arg[iarg+1], "nohistory") == 0){ + if (twist != TWIST_NOHISTORY) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be the same for all types"); } - else if (strstr(arg[iarg+1], "sds") != NULL){ - if (twist != TWIST_SDS) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be consistent"); + else if (strcmp(arg[iarg+1], "sds") == 0){ + if (!twist_set) twist = TWIST_SDS; + else if (twist != TWIST_SDS) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be the same for all types"); twist_history = 1; } else{ @@ -1727,6 +1548,15 @@ void PairGranular::coeff(int narg, char **arg) else error->all(FLERR, "Illegal pair coeff command"); } + //It is an error not to specify normal or tangential model + if (!normal_set) error->all(FLERR, "Illegal pair coeff command, must specify normal contact model"); + if (!tangential_set) error->all(FLERR, "Illegal pair coeff command, must specify tangential contact model"); + + //If unspecified, set damping to VISCOELASTIC, twist/roll to NONE (cannot be changed by subsequent pair_coeff commands) + if (!damping_set) damping = VISCOELASTIC; + if (!roll_set) roll = ROLL_NONE; + if (!twist_set) twist = TWIST_NONE; + int count = 0; double damp; if (damping == TSUJI){ @@ -1740,37 +1570,28 @@ void PairGranular::coeff(int narg, char **arg) for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { - if (normal_set){ - normal_coeffs[i][j][0] = normal_coeffs_local[0]; - normal_coeffs[i][j][1] = damp; - if (normal != HERTZ && normal != HOOKE) normal_coeffs[i][j][2] = normal_coeffs_local[2]; - if ((normal == JKR) || (normal == DMT)) - normal_coeffs[i][j][3] = normal_coeffs_local[3]; - } - if (tangential_set){ + normal_coeffs[i][j][0] = normal_coeffs_local[0]; + normal_coeffs[i][j][1] = damp; + if (normal != HERTZ && normal != HOOKE) normal_coeffs[i][j][2] = normal_coeffs_local[2]; + if ((normal == JKR) || (normal == DMT)) + normal_coeffs[i][j][3] = normal_coeffs_local[3]; + + for (int k = 0; k < 3; k++) + tangential_coeffs[i][j][k] = tangential_coeffs_local[k]; + + if (roll != ROLL_NONE) for (int k = 0; k < 3; k++) - tangential_coeffs[i][j][k] = tangential_coeffs_local[k]; - } - if (roll_set){ - if (roll != ROLL_NONE) - for (int k = 0; k < 3; k++) - roll_coeffs[i][j][k] = roll_coeffs_local[k]; - } - if (twist_set){ - if (twist != TWIST_NONE && twist != TWIST_MARSHALL) - for (int k = 0; k < 3; k++) - twist_coeffs[i][j][k] = twist_coeffs_local[k]; - } + roll_coeffs[i][j][k] = roll_coeffs_local[k]; + + if (twist != TWIST_NONE && twist != TWIST_MARSHALL) + for (int k = 0; k < 3; k++) + twist_coeffs[i][j][k] = twist_coeffs_local[k]; + setflag[i][j] = 1; count++; } } - delete[] normal_coeffs_local; - delete[] tangential_coeffs_local; - delete[] roll_coeffs_local; - delete[] twist_coeffs_local; - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } @@ -1959,20 +1780,25 @@ double PairGranular::init_one(int i, int j) // if there is no current information about radius/cutoff of type i and j). // we assign cutoff = max(cut[i][j]) for i,j such that cut[i][j] > 0.0. - if (((maxrad_dynamic[i] > 0.0) && (maxrad_dynamic[j] > 0.0)) || - ((maxrad_dynamic[i] > 0.0) && (maxrad_frozen[j] > 0.0)) || - ((maxrad_frozen[i] > 0.0) && (maxrad_dynamic[j] > 0.0))) { // radius info about both i and j exist - cutoff = maxrad_dynamic[i]+maxrad_dynamic[j]; - cutoff = MAX(cutoff,maxrad_frozen[i]+maxrad_dynamic[j]); - cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]); - } - else { // radius info about either i or j does not exist (i.e. not present and not about to get poured; set to largest value to not interfere with neighbor list) - double cutmax = 0.0; - for (int k = 1; k <= atom->ntypes; k++) { - cutmax = MAX(cutmax,2.0*maxrad_dynamic[k]); - cutmax = MAX(cutmax,2.0*maxrad_frozen[k]); + if (cutoff_global < 0){ + if (((maxrad_dynamic[i] > 0.0) && (maxrad_dynamic[j] > 0.0)) || + ((maxrad_dynamic[i] > 0.0) && (maxrad_frozen[j] > 0.0)) || + ((maxrad_frozen[i] > 0.0) && (maxrad_dynamic[j] > 0.0))) { // radius info about both i and j exist + cutoff = maxrad_dynamic[i]+maxrad_dynamic[j]; + cutoff = MAX(cutoff,maxrad_frozen[i]+maxrad_dynamic[j]); + cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]); } - cutoff = cutmax; + else { // radius info about either i or j does not exist (i.e. not present and not about to get poured; set to largest value to not interfere with neighbor list) + double cutmax = 0.0; + for (int k = 1; k <= atom->ntypes; k++) { + cutmax = MAX(cutmax,2.0*maxrad_dynamic[k]); + cutmax = MAX(cutmax,2.0*maxrad_frozen[k]); + } + cutoff = cutmax; + } + } + else{ + cutoff = cutoff_global; } return cutoff; } diff --git a/src/GRANULAR/pair_granular.h b/src/GRANULAR/pair_granular.h index 897316c907..f3a9d4dcbe 100644 --- a/src/GRANULAR/pair_granular.h +++ b/src/GRANULAR/pair_granular.h @@ -86,19 +86,18 @@ private: //Indices of history entries int tangential_history_index, roll_history_index, twist_history_index; - //Coefficients declared in pair style command, used as default unless - // overwritten in pair coeff command - double *normal_coeffs_global; - double *tangential_coeffs_global; - double *roll_coeffs_global; - double *twist_coeffs_global; + //Flags for whether model choices have been set + int normal_set, tangential_set, damping_set, roll_set, twist_set; - //Per-type coefficients declared in pair coeff command + //Per-type coefficients, set in pair coeff command double ***normal_coeffs; double ***tangential_coeffs; double ***roll_coeffs; double ***twist_coeffs; + //Optional user-specified global cutoff + double global_cutoff; + double mix_stiffnessE(double Eii, double Ejj, double Gii, double Gjj); double mix_stiffnessG(double Eii, double Ejj, double Gii, double Gjj); double mix_geom(double valii, double valjj); diff --git a/src/GRANULAR/pair_granular_multi.cpp b/src/GRANULAR/pair_granular_multi.cpp index 11381444a2..5fee363872 100644 --- a/src/GRANULAR/pair_granular_multi.cpp +++ b/src/GRANULAR/pair_granular_multi.cpp @@ -654,217 +654,15 @@ void PairGranularMulti::allocate() void PairGranularMulti::settings(int narg, char **arg) { - if (narg < 5) error->all(FLERR,"Illegal pair_style command"); + if (narg == 1){ + cutoff_global = force->numeric(FLERR,arg[0]) + } + else{ + cutoff_global = -1; //Will be set based on particle sizes, model choice + } - int iarg = 0; - - //Some defaults - normal_global = HERTZ; - damping_global = VISCOELASTIC; - tangential_global = TANGENTIAL_MINDLIN; - roll_global = ROLL_NONE; - twist_global = TWIST_NONE; - - tangential_history = 1; + tangential_history = 0; roll_history = twist_history = 0; - - int normal_set, tangential_set; - normal_set = tangential_set = 0; - - while (iarg < narg){ - if (strcmp(arg[iarg], "hooke") == 0){ - if (iarg + 2 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for Hooke option"); - normal_global = HOOKE; - normal_set = 1; - memory->create(normal_coeffs_global, 2, "pair:normal_coeffs_global"); - normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1]); //kn - normal_coeffs_global[1] = force->numeric(FLERR,arg[iarg+2]); //damping - iarg += 3; - } - else if (strcmp(arg[iarg], "hertz") == 0){ - int num_coeffs = 2; - if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for Hertz option"); - normal_global = HERTZ; - normal_set = 1; - memory->create(normal_coeffs_global, num_coeffs, "pair:normal_coeffs_global"); - normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1]); //kn - normal_coeffs_global[1] = force->numeric(FLERR,arg[iarg+2]); //damping - iarg += num_coeffs+1; - } - else if (strcmp(arg[iarg], "hertz/material") == 0){ - int num_coeffs = 3; - if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for Hertz option"); - normal_global = HERTZ_MATERIAL; - normal_set = 1; - memory->create(normal_coeffs_global, num_coeffs, "pair:normal_coeffs_global"); - normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //E (Young's modulus) - normal_coeffs_global[1] = force->numeric(FLERR,arg[iarg+2]); //damping - normal_coeffs_global[2] = force->numeric(FLERR,arg[iarg+3]); //G (shear modulus) - iarg += num_coeffs+1; - } - else if (strcmp(arg[iarg], "dmt") == 0){ - if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for Hertz option"); - normal_global = DMT; - normal_set = 1; - memory->create(normal_coeffs_global, 4, "pair:normal_coeffs_global"); - normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //4/3 E - normal_coeffs_global[1] = force->numeric(FLERR,arg[iarg+2]); //damping - normal_coeffs_global[2] = force->numeric(FLERR,arg[iarg+3]); //G - normal_coeffs_global[3] = force->numeric(FLERR,arg[iarg+3]); //cohesion - iarg += 5; - } - else if (strcmp(arg[iarg], "jkr") == 0){ - if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for JKR option"); - beyond_contact = 1; - normal_global = JKR; - normal_set = 1; - memory->create(normal_coeffs_global, 4, "pair:normal_coeffs_global"); - normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1]); //E - normal_coeffs_global[1] = force->numeric(FLERR,arg[iarg+2]); //damping - normal_coeffs_global[2] = force->numeric(FLERR,arg[iarg+3]); //G - normal_coeffs_global[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion - iarg += 5; - } - else if (strcmp(arg[iarg], "damp_velocity") == 0){ - damping_global = VELOCITY; - iarg += 1; - } - else if (strcmp(arg[iarg], "damp_viscoelastic") == 0){ - damping_global = VISCOELASTIC; - iarg += 1; - } - else if (strcmp(arg[iarg], "damp_tsuji") == 0){ - damping_global = TSUJI; - iarg += 1; - } - else if (strstr(arg[iarg], "tangential") != NULL){ - if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for tangential model"); - if (strstr(arg[iarg+1], "nohistory") != NULL){ - tangential_global = TANGENTIAL_NOHISTORY; - tangential_set = 1; - } - else if (strstr(arg[iarg+1], "mindlin") != NULL){ - tangential_global = TANGENTIAL_MINDLIN; - tangential_history = 1; - tangential_set = 1; - } - else{ - error->all(FLERR, "Illegal pair_style command, unrecognized sliding friction model"); - } - memory->create(tangential_coeffs_global, 3, "pair:tangential_coeffs_global"); - tangential_coeffs_global[0] = force->numeric(FLERR,arg[iarg+2]); //kt - tangential_coeffs_global[1] = force->numeric(FLERR,arg[iarg+3]); //gammat - tangential_coeffs_global[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. - iarg += 5; - } - else if (strstr(arg[iarg], "roll") != NULL){ - if (strstr(arg[iarg+1], "none") != NULL){ - roll_global = ROLL_NONE; - iarg += 2; - } - else{ - if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for rolling model"); - if (strstr(arg[iarg+1], "nohistory") != NULL){ - roll_global = ROLL_NOHISTORY; - } - else if (strstr(arg[iarg+1], "sds") != NULL){ - roll_global = ROLL_SDS; - roll_history = 1; - } - else{ - error->all(FLERR, "Illegal pair_style command, unrecognized rolling friction model"); - } - memory->create(roll_coeffs_global, 3, "pair:roll_coeffs_global"); - roll_coeffs_global[0] = force->numeric(FLERR,arg[iarg+2]); //kR - roll_coeffs_global[1] = force->numeric(FLERR,arg[iarg+3]); //gammaR - roll_coeffs_global[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. - iarg += 5; - } - } - else if (strstr(arg[iarg], "twist") != NULL){ - if (strstr(arg[iarg+1], "none") != NULL){ - twist_global = TWIST_NONE; - iarg += 2; - } - else if (strstr(arg[iarg+1], "marshall") != NULL){ - twist_global = TWIST_MARSHALL; - twist_history = 1; - iarg += 2; - } - else{ - if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for twist model"); - memory->create(twist_coeffs_global, 3, "pair:twist_coeffs_global"); //To be filled later - if (strstr(arg[iarg+1], "nohistory") != NULL){ - twist_global = TWIST_NOHISTORY; - } - else if (strstr(arg[iarg+1], "sds") != NULL){ - twist_global = TWIST_SDS; - twist_history = 1; - } - else{ - error->all(FLERR, "Illegal pair_style command, unrecognized twisting friction model"); - } - memory->create(twist_coeffs_global, 3, "pair:twist_coeffs_global"); - twist_coeffs_global[0] = force->numeric(FLERR,arg[iarg+2]); //ktwist - twist_coeffs_global[1] = force->numeric(FLERR,arg[iarg+3]); //gammatwist - twist_coeffs_global[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. - iarg += 5; - } - } - else error->all(FLERR, "Illegal pair_style granular command"); - } - - //Set all i-i entries, which may be replaced by pair coeff commands - //It may also make sense to consider removing all of the above, and only - // having the option for pair_coeff to set the parameters, similar to most LAMMPS pair styles - // The reason for the current setup is to remain true to existing pair gran/hooke etc. syntax, - // where coeffs are set in the pair_style command, and a pair_coeff * * command is issued. - - //Other option is to have two pair styles, e.g. pair granular and pair granular/multi, - // where granular/multi allows per-type coefficients, pair granular does not (this would also - // allow minor speed-up by templating pair granular) - allocate(); - double damp; - for (int i = 1; i <= atom->ntypes; i++){ - normal[i][i] = normal_global; - damping[i][i] = damping_global; - tangential[i][i] = tangential_global; - roll[i][i] = roll_global; - twist[i][i] = twist_global; - - if (normal_set){ - if (damping_global == TSUJI){ - double cor = normal_coeffs_global[1]; - damp = 1.2728-4.2783*cor+11.087*pow(cor,2)-22.348*pow(cor,3)+ - 27.467*pow(cor,4)-18.022*pow(cor,5)+ - 4.8218*pow(cor,6); - } - else damp = normal_coeffs_global[1]; - normal_coeffs[i][i][0] = normal_coeffs_global[0]; - normal_coeffs[i][i][1] = damp; - if (normal[i][i] != HOOKE && normal[i][i] != HERTZ){ - normal_coeffs[i][i][2] = normal_coeffs_global[2]; - } - if ((normal_global == JKR) || (normal_global == DMT)) - normal_coeffs[i][i][3] = normal_coeffs_global[3]; - } - if(tangential_set){ - tangential[i][i] = tangential_global; - for (int k = 0; k < 3; k++) - tangential_coeffs[i][i][k] = tangential_coeffs_global[k]; - } - roll[i][i] = roll_global; - if (roll_global != ROLL_NONE) - for (int k = 0; k < 3; k++) - roll_coeffs[i][i][k] = roll_coeffs_global[k]; - - twist[i][i] = twist_global; - if (twist_global != TWIST_NONE && twist_global != TWIST_MARSHALL) - for (int k = 0; k < 3; k++) - twist_coeffs[i][i][k] = twist_coeffs_global[k]; - - if (normal_set && tangential_set) setflag[i][i] = 1; - } } /* ---------------------------------------------------------------------- @@ -874,10 +672,6 @@ void PairGranularMulti::settings(int narg, char **arg) void PairGranularMulti::coeff(int narg, char **arg) { int normal_local, damping_local, tangential_local, roll_local, twist_local; - double *normal_coeffs_local; - double *tangential_coeffs_local; - double *roll_coeffs_local; - double *twist_coeffs_local; normal_coeffs_local = new double[4]; tangential_coeffs_local = new double[4]; @@ -893,8 +687,10 @@ void PairGranularMulti::coeff(int narg, char **arg) force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi); force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi); - normal_local = tangential_local = roll_local = twist_local = -1; - damping_local = -1; + //Defaults + normal_local = tangential_local = -1; + roll_local = twist_local = 0; + damping_local = VISCOELASTIC; int iarg = 2; while (iarg < narg){ @@ -941,24 +737,27 @@ void PairGranularMulti::coeff(int narg, char **arg) normal_coeffs_local[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion iarg += 5; } - else if (strcmp(arg[iarg], "damp_velocity") == 0){ - damping_local = VELOCITY; - iarg += 1; + else if (strcmp(arg[iarg], "damp") == 0){ + if (iarg+1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters provided for damping model"); + if (strcmp(arg[iarg+1]), "velocity") == 0){ + damping_local = VELOCITY; + iarg += 1; + } + else if (strcmp(arg[iarg+1], "viscoelastic") == 0){ + damping_local = VISCOELASTIC; + iarg += 1; + } + else if (strcmp(arg[iarg], "tsuji") == 0){ + damping_local = TSUJI; + iarg += 1; + } } - else if (strcmp(arg[iarg], "damp_viscoelastic") == 0){ - damping_local = VISCOELASTIC; - iarg += 1; - } - else if (strcmp(arg[iarg], "damp_tsuji") == 0){ - damping_local = TSUJI; - iarg += 1; - } - else if (strstr(arg[iarg], "tangential") != NULL){ + else if (strcmp(arg[iarg], "tangential") == 0){ if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model"); - if (strstr(arg[iarg+1], "nohistory") != NULL){ + if (strcmp(arg[iarg+1], "nohistory") == 0){ tangential_local = TANGENTIAL_NOHISTORY; } - else if (strstr(arg[iarg+1], "mindlin") != NULL){ + else if (strcmp(arg[iarg+1], "mindlin") == 0){ tangential_local = TANGENTIAL_MINDLIN; tangential_history = 1; } @@ -970,18 +769,18 @@ void PairGranularMulti::coeff(int narg, char **arg) tangential_coeffs_local[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. iarg += 5; } - else if (strstr(arg[iarg], "rolling") != NULL){ + else if (strcmp(arg[iarg], "rolling") == 0){ if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); - if (strstr(arg[iarg+1], "none") != NULL){ + if (strcmp(arg[iarg+1], "none") == 0){ roll_local = ROLL_NONE; iarg += 2; } else{ if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for rolling model"); - if (strstr(arg[iarg+1], "nohistory") != NULL){ + if (strcmp(arg[iarg+1], "nohistory") == 0){ roll_local = ROLL_NOHISTORY; } - else if (strstr(arg[iarg+1], "sds") != NULL){ + else if (strcmp(arg[iarg+1], "sds") == 0){ roll_local = ROLL_SDS; roll_history = 1; } @@ -994,23 +793,23 @@ void PairGranularMulti::coeff(int narg, char **arg) iarg += 5; } } - else if (strstr(arg[iarg], "twist") != NULL){ + else if (strcmp(arg[iarg], "twisting") == 0){ if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); - if (strstr(arg[iarg+1], "none") != NULL){ + if (strcmp(arg[iarg+1], "none") == 0){ twist_local = TWIST_NONE; iarg += 2; } - else if (strstr(arg[iarg+1], "marshall") != NULL){ + else if (strcmp(arg[iarg+1], "marshall") == 0){ twist_local = TWIST_MARSHALL; twist_history = 1; iarg += 2; } else{ if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for twist model"); - if (strstr(arg[iarg+1], "nohistory") != NULL){ + if (strcmp(arg[iarg+1], "nohistory") == 0){ twist_local = TWIST_NOHISTORY; } - else if (strstr(arg[iarg+1], "sds") != NULL){ + else if (strcmp(arg[iarg+1], "sds") == 0){ twist_local = TWIST_SDS; twist_history = 1; } @@ -1026,66 +825,49 @@ void PairGranularMulti::coeff(int narg, char **arg) else error->all(FLERR, "Illegal pair coeff command"); } + //It is an error not to specify normal or tangential model + if ((normal_set < 0) || (tangential_set < 0)) error->all(FLERR, "Illegal pair coeff command, must specify normal contact model");)) + int count = 0; double damp; - if (damping_local >= 0){ - if (normal_local == -1) - error->all(FLERR, "Illegal pair_coeff command, must specify normal model when setting damping model"); - if (damping_local == TSUJI){ - double cor; - cor = normal_coeffs_local[1]; - damp = 1.2728-4.2783*cor+11.087*pow(cor,2)-22.348*pow(cor,3)+ - 27.467*pow(cor,4)-18.022*pow(cor,5)+ - 4.8218*pow(cor,6); - } - else damp = normal_coeffs_local[1]; + if (damping_local == TSUJI){ + double cor; + cor = normal_coeffs_local[1]; + damp = 1.2728-4.2783*cor+11.087*pow(cor,2)-22.348*pow(cor,3)+ + 27.467*pow(cor,4)-18.022*pow(cor,5)+ + 4.8218*pow(cor,6); } + else damp = normal_coeffs_local[1]; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { - if (normal_local >= 0){ - normal[i][j] = normal_local; - normal_coeffs[i][j][0] = normal_coeffs_local[0]; - if (damping_local == -1){ - damp = normal_coeffs_global[1]; - } - normal_coeffs[i][j][1] = damp; - if (normal_local != HERTZ && normal_local != HOOKE) normal_coeffs[i][j][2] = normal_coeffs_local[2]; - if ((normal_local == JKR) || (normal_local == DMT)) - normal_coeffs[i][j][3] = normal_coeffs_local[3]; - } - if (damping_local >= 0){ - damping[i][j] = damping_local; - } - if (tangential_local >= 0){ - tangential[i][j] = tangential_local; + normal[i][j] = normal_local; + normal_coeffs[i][j][0] = normal_coeffs_local[0]; + normal_coeffs[i][j][1] = damp; + if (normal_local != HERTZ && normal_local != HOOKE) normal_coeffs[i][j][2] = normal_coeffs_local[2]; + if ((normal_local == JKR) || (normal_local == DMT)) + normal_coeffs[i][j][3] = normal_coeffs_local[3]; + + damping[i][j] = damping_local; + + tangential[i][j] = tangential_local; for (int k = 0; k < 3; k++) tangential_coeffs[i][j][k] = tangential_coeffs_local[k]; - } - if (roll_local >= 0){ - roll[i][j] = roll_local; - if (roll_local != ROLL_NONE) - for (int k = 0; k < 3; k++) - roll_coeffs[i][j][k] = roll_coeffs_local[k]; - } - if (twist_local >= 0){ - twist[i][j] = twist_local; + + roll[i][j] = roll_local; + if (roll_local != ROLL_NONE) + for (int k = 0; k < 3; k++) + roll_coeffs[i][j][k] = roll_coeffs_local[k]; + + twist[i][j] = twist_local; if (twist_local != TWIST_NONE && twist_local != TWIST_MARSHALL) for (int k = 0; k < 3; k++) twist_coeffs[i][j][k] = twist_coeffs_local[k]; - } - - if (normal_local >= 0 && tangential_local >= 0) setflag[i][j] = 1; + setflag[i][j] = 1; count++; } } - - delete[] normal_coeffs_local; - delete[] tangential_coeffs_local; - delete[] roll_coeffs_local; - delete[] twist_coeffs_local; - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } @@ -1286,20 +1068,25 @@ double PairGranularMulti::init_one(int i, int j) // if there is no current information about radius/cutoff of type i and j). // we assign cutoff = max(cut[i][j]) for i,j such that cut[i][j] > 0.0. - if (((maxrad_dynamic[i] > 0.0) && (maxrad_dynamic[j] > 0.0)) || - ((maxrad_dynamic[i] > 0.0) && (maxrad_frozen[j] > 0.0)) || - ((maxrad_frozen[i] > 0.0) && (maxrad_dynamic[j] > 0.0))) { // radius info about both i and j exist - cutoff = maxrad_dynamic[i]+maxrad_dynamic[j]; - cutoff = MAX(cutoff,maxrad_frozen[i]+maxrad_dynamic[j]); - cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]); - } - else { // radius info about either i or j does not exist (i.e. not present and not about to get poured; set to largest value to not interfere with neighbor list) - double cutmax = 0.0; - for (int k = 1; k <= atom->ntypes; k++) { - cutmax = MAX(cutmax,2.0*maxrad_dynamic[k]); - cutmax = MAX(cutmax,2.0*maxrad_frozen[k]); + if (cutoff_global < 0){ + if (((maxrad_dynamic[i] > 0.0) && (maxrad_dynamic[j] > 0.0)) || + ((maxrad_dynamic[i] > 0.0) && (maxrad_frozen[j] > 0.0)) || + ((maxrad_frozen[i] > 0.0) && (maxrad_dynamic[j] > 0.0))) { // radius info about both i and j exist + cutoff = maxrad_dynamic[i]+maxrad_dynamic[j]; + cutoff = MAX(cutoff,maxrad_frozen[i]+maxrad_dynamic[j]); + cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]); } - cutoff = cutmax; + else { // radius info about either i or j does not exist (i.e. not present and not about to get poured; set to largest value to not interfere with neighbor list) + double cutmax = 0.0; + for (int k = 1; k <= atom->ntypes; k++) { + cutmax = MAX(cutmax,2.0*maxrad_dynamic[k]); + cutmax = MAX(cutmax,2.0*maxrad_frozen[k]); + } + cutoff = cutmax; + } + } + else{ + cutoff = cutoff_global; } return cutoff; } diff --git a/src/GRANULAR/pair_granular_multi.h b/src/GRANULAR/pair_granular_multi.h index e853e564df..95beb950f4 100644 --- a/src/GRANULAR/pair_granular_multi.h +++ b/src/GRANULAR/pair_granular_multi.h @@ -65,27 +65,26 @@ public: private: int size_history; - //Per-type models + //Models int **normal, **damping, **tangential, **roll, **twist; - int normal_global, damping_global; - int tangential_global, roll_global, twist_global; - + //History flags int tangential_history, roll_history, twist_history; + + //Indices of history entries int tangential_history_index; int roll_history_index; int twist_history_index; - double *normal_coeffs_global; - double *tangential_coeffs_global; - double *roll_coeffs_global; - double *twist_coeffs_global; - + //Per-type coefficients, set in pair coeff command double ***normal_coeffs; double ***tangential_coeffs; double ***roll_coeffs; double ***twist_coeffs; + //Optional user-specified global cutoff + double cutoff_global; + double mix_stiffnessE(double Eii, double Ejj, double Gii, double Gjj); double mix_stiffnessG(double Eii, double Ejj, double Gii, double Gjj); double mix_geom(double valii, double valjj);