From 038f4a52105aa93bbd835310d65be697812109f8 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Wed, 14 Sep 2022 21:40:00 -0600 Subject: [PATCH] Adding variable temperature to fix gran/wall, misc fixes/updates --- src/GRANULAR/contact.cpp | 62 ++++++++-------------- src/GRANULAR/contact.h | 12 ++--- src/GRANULAR/contact_normal_models.cpp | 2 +- src/GRANULAR/contact_sub_models.cpp | 12 +++-- src/GRANULAR/contact_tangential_models.cpp | 2 +- src/GRANULAR/contact_twisting_models.cpp | 16 ++++-- src/GRANULAR/contact_twisting_models.h | 3 ++ src/GRANULAR/fix_wall_gran.cpp | 29 +++++++--- src/GRANULAR/fix_wall_gran.h | 3 ++ src/GRANULAR/fix_wall_gran_region.cpp | 5 +- src/GRANULAR/pair_granular.cpp | 33 ++++++------ 11 files changed, 98 insertions(+), 81 deletions(-) diff --git a/src/GRANULAR/contact.cpp b/src/GRANULAR/contact.cpp index a1a35afe47..f72245def3 100644 --- a/src/GRANULAR/contact.cpp +++ b/src/GRANULAR/contact.cpp @@ -10,12 +10,16 @@ the GNU General Public License. See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- +------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- This class contains a series of tools for DEM contacts Multiple models can be defined and used to calculate forces and torques based on contact geometry -*/ + + Contributing authors: + Dan Bolintineanu (SNL), Joel Clemmer (SNL) +----------------------------------------------------------------------- */ #include "contact.h" #include "contact_sub_models.h" @@ -42,7 +46,7 @@ ContactModel::ContactModel(LAMMPS *lmp) : Pointers(lmp) beyond_contact = 0; nondefault_history_transfer = 0; - wall_type = NONE; + contact_type = PAIR; reset_contact(); @@ -175,17 +179,19 @@ int ContactModel::init_classic_model(char **arg, int iarg, int narg) if (strcmp(arg[iarg],"hooke") == 0) { init_model("hooke", NORMAL); init_model("linear_nohistory", TANGENTIAL); + init_model("mass_velocity", DAMPING); } else if (strcmp(arg[iarg],"hooke/history") == 0) { init_model("hooke", NORMAL); init_model("linear_history", TANGENTIAL); + init_model("mass_velocity", DAMPING); } else if (strcmp(arg[iarg],"hertz/history") == 0) { // convert Kn and Kt from pressure units to force/distance^2 if Hertzian kn /= force->nktv2p; kt /= force->nktv2p; init_model("hertz", NORMAL); - init_model("mindlin", TANGENTIAL); // Dan is this right? + init_model("mindlin", TANGENTIAL); + init_model("viscoelastic", DAMPING); } else error->all(FLERR,"Invalid classic gran model"); - init_model("mass_velocity", DAMPING); // ensure additional models are undefined init_model("none", ROLLING); @@ -322,25 +328,15 @@ void ContactModel::read_restart(FILE *fp) /* ---------------------------------------------------------------------- */ -void ContactModel::reset_contact() -{ - prep_flag = check_flag = 0; - touch = false; -} - -/* ---------------------------------------------------------------------- */ - bool ContactModel::check_contact(double rtemp) { - check_flag = 1; - - if (wall_type == RWALL) { + if (contact_type == WALL) { // Used by fix_wall_gran.cpp rsq = lensq3(dx); radsum = radi; if (rtemp == 0) Reff = radi; else Reff = radi * rtemp/(radi + rtemp); - } else if (wall_type == RDUPLICATE) { + } else if (contact_type == WALLREGION) { // Used by fix_wall_gran_region.cpp rsq = rtemp * rtemp; radsum = radi + radi; @@ -361,11 +357,6 @@ bool ContactModel::check_contact(double rtemp) void ContactModel::prep_contact() { - prep_flag = 1; - // If it hasn't already been done, test if the contact exists - if (check_flag != 1) touch = check_contact(); - if (!touch) return; - double temp[3]; // Standard geometric quantities @@ -415,16 +406,8 @@ void ContactModel::prep_contact() void ContactModel::calculate_forces() { - // If it hasn't already been done, run prep calculations - if (prep_flag != 1) prep_contact(); - if (!touch) { - forces[0] = forces[1] = forces[2] = 0.0; - return; - } + // calculate forces/torques - //********************************************** - // calculate forces - //********************************************** forces[0] = 0.0; double Fne, Fdamp; area = normal_model->calculate_area(); @@ -442,35 +425,36 @@ void ContactModel::calculate_forces() if (rolling_defined) rolling_model->calculate_forces(); if (twisting_defined) twisting_model->calculate_forces(); - //********************************************** // sum contributions - //********************************************** scale3(Fntot, nx, forces); add3(forces, fs, forces); - //May need to rethink this for use with walls (and eventually tris).. + //May need to rethink eventually tris.. cross3(nx, fs, torquesi); - copy3(torquesi, torquesj); double dist_to_contact = radi - 0.5 * delta; scale3(dist_to_contact, torquesi); - dist_to_contact = radj - 0.5 * delta; - scale3(dist_to_contact, torquesj); + + if (contact_type == PAIR) { + copy3(torquesi, torquesj); + dist_to_contact = radj - 0.5 * delta; + scale3(dist_to_contact, torquesj); + } double torroll[3]; if (rolling_defined) { cross3(nx, fr, torroll); scale3(Reff, torroll); add3(torquesi, torroll, torquesi); - sub3(torquesj, torroll, torquesj); + if (contact_type == PAIR) sub3(torquesj, torroll, torquesj); } double tortwist[3]; if (twisting_defined) { scale3(magtortwist, nx, tortwist); add3(torquesi, tortwist, torquesi); - sub3(torquesj, tortwist, torquesj); + if (contact_type == PAIR) sub3(torquesj, tortwist, torquesj); } if (heat_defined) { diff --git a/src/GRANULAR/contact.h b/src/GRANULAR/contact.h index e0600146c3..4dc88132c0 100644 --- a/src/GRANULAR/contact.h +++ b/src/GRANULAR/contact.h @@ -31,10 +31,10 @@ enum ModelType { HEAT = 5 }; // Relative order matters since some derive coeffs from others -enum WallType { - NONE = 0, - RWALL = 1, - RDUPLICATE = 2 +enum ContactType { + PAIR = 0, + WALL = 1, + WALLREGION = 2 }; // forward declaration @@ -52,7 +52,6 @@ class ContactModel : protected Pointers { ~ContactModel(); void init(); bool check_contact(double = 0); - void reset_contact(); void prep_contact(); void calculate_forces(); double pulloff_distance(double, double); @@ -75,7 +74,7 @@ class ContactModel : protected Pointers { // Extra options int beyond_contact, limit_damping, history_update; - WallType wall_type; + ContactType contact_type; // History variables int size_history, nondefault_history_transfer; @@ -98,7 +97,6 @@ class ContactModel : protected Pointers { protected: int rolling_defined, twisting_defined, heat_defined; // Used to quickly skip undefined submodels - int prep_flag, check_flag; }; } // namespace Contact diff --git a/src/GRANULAR/contact_normal_models.cpp b/src/GRANULAR/contact_normal_models.cpp index 2f59311edf..9e0a377eb9 100644 --- a/src/GRANULAR/contact_normal_models.cpp +++ b/src/GRANULAR/contact_normal_models.cpp @@ -261,7 +261,7 @@ void NormalJKR::coeffs_to_local() poiss = coeffs[2]; cohesion = coeffs[3]; k = FOURTHIRDS * Emod; - Escaled = mix_stiffnessE(Emod, Emod, poiss, poiss); //Dan, not sure why these coefficients are mixed in the regular pair style + Escaled = mix_stiffnessE(Emod, Emod, poiss, poiss); if (Emod < 0.0 || damp < 0.0) error->all(FLERR, "Illegal JKR normal model"); } diff --git a/src/GRANULAR/contact_sub_models.cpp b/src/GRANULAR/contact_sub_models.cpp index 6ec7bbd017..57834fc040 100644 --- a/src/GRANULAR/contact_sub_models.cpp +++ b/src/GRANULAR/contact_sub_models.cpp @@ -10,12 +10,16 @@ the GNU General Public License. See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- +-------------------------------------------------------------------------*/ - This class contains a series of tools for DEM contacts - Multiple models can be defined and used to calculate forces +/* ---------------------------------------------------------------------- + This class contains a framework for normal, damping, tangential, + rolling, twisting, and heat models used to calculate forces and torques based on contact geometry -*/ + + Contributing authors: + Dan Bolintineanu (SNL), Joel Clemmer (SNL) +----------------------------------------------------------------------- */ #include "contact_sub_models.h" #include "error.h" diff --git a/src/GRANULAR/contact_tangential_models.cpp b/src/GRANULAR/contact_tangential_models.cpp index 437615441c..46b370d15f 100644 --- a/src/GRANULAR/contact_tangential_models.cpp +++ b/src/GRANULAR/contact_tangential_models.cpp @@ -259,7 +259,7 @@ void TangentialMindlin::calculate_forces() temp_dbl = -damp; scale3(temp_dbl, contact->vtr, contact->fs); - if (! mindlin_force) { + if (!mindlin_force) { scale3(k_scaled, history, temp_array); add3(contact->fs, temp_array, contact->fs); } diff --git a/src/GRANULAR/contact_twisting_models.cpp b/src/GRANULAR/contact_twisting_models.cpp index 8f3c5e571e..52b37d46e1 100644 --- a/src/GRANULAR/contact_twisting_models.cpp +++ b/src/GRANULAR/contact_twisting_models.cpp @@ -46,15 +46,25 @@ TwistingMarshall::TwistingMarshall(LAMMPS *lmp) : TwistingModel(lmp) /* ---------------------------------------------------------------------- */ + +void TwistingMarshall::init() +{ + k_tang = contact->tangential_model->k; + damp_tang = contact->tangential_model->damp; + mu_tang = contact->tangential_model->mu; +} + +/* ---------------------------------------------------------------------- */ + void TwistingMarshall::calculate_forces() { double signtwist, Mtcrit; // Calculate twist coefficients from tangential model & contact geometry // eq 32 of Marshall paper - double k = 0.5 * contact->tangential_model->k * contact->area * contact->area; - double damp = 0.5 * contact->tangential_model->damp * contact->area * contact->area; - double mu = TWOTHIRDS * contact->area * contact->tangential_model->mu; + double k = 0.5 * k_tang * contact->area * contact->area; + double damp = 0.5 * damp_tang * contact->area * contact->area; + double mu = TWOTHIRDS * mu_tang * contact->area; if (contact->history_update) { contact->history[history_index] += contact->magtwist * contact->dt; diff --git a/src/GRANULAR/contact_twisting_models.h b/src/GRANULAR/contact_twisting_models.h index 15124ca378..724448892a 100644 --- a/src/GRANULAR/contact_twisting_models.h +++ b/src/GRANULAR/contact_twisting_models.h @@ -42,6 +42,9 @@ class TwistingMarshall : public TwistingModel { public: TwistingMarshall(class LAMMPS *); void calculate_forces(); + void init(); + protected: + double k_tang, damp_tang, mu_tang; }; /* ---------------------------------------------------------------------- */ diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index 95cd4fc6e0..aca9341a58 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -14,7 +14,7 @@ /* ---------------------------------------------------------------------- Contributing authors: Leo Silbert (SNL), Gary Grest (SNL), - Dan Bolintineanu (SNL) + Dan Bolintineanu (SNL), Joel Clemmer (SNL) ------------------------------------------------------------------------- */ #include "fix_wall_gran.h" @@ -69,7 +69,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : // set interaction style // disable bonded/history option for now model = new ContactModel(lmp); - model->wall_type = RWALL; + model->contact_type = WALL; if (strcmp(arg[3],"granular") == 0) classic_flag = 0; else classic_flag = 1; @@ -171,6 +171,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : // wallstyle args idregion = nullptr; + tstr = nullptr; if (iarg >= narg) error->all(FLERR, "Illegal fix wall/gran command"); @@ -211,7 +212,11 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : iarg += 2; } else if (strcmp(arg[iarg],"temperature") == 0) { if (narg < iarg+2) error->all(FLERR,"Illegal fix wall/gran command"); - Twall = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (utils::strmatch(arg[iarg+1], "^v_")) { + tstr = utils::strdup(arg[3] + 2); + } else { + Twall = utils::numeric(FLERR,arg[iarg+1],false,lmp); + } Twall_defined = 1; iarg += 2; } else wallstyle = NOSTYLE; @@ -318,10 +323,10 @@ FixWallGran::~FixWallGran() atom->delete_callback(id,Atom::GROW); atom->delete_callback(id,Atom::RESTART); - delete model; - // delete local storage + delete model; + delete [] tstr; delete [] idregion; memory->destroy(history_one); memory->destroy(mass_rigid); @@ -365,6 +370,13 @@ void FixWallGran::init() model->sub_models[i]->history_index = next_index; next_index += model->sub_models[i]->size_history; } + + if (tstr) { + tvar = input->variable->find(tstr); + if (tvar < 0) error->all(FLERR, "Variable {} for fix wall/gran does not exist", tstr); + if (! input->variable->equalstyle(tvar)) + error->all(FLERR, "Variable {} for fix wall/gran must be an equal style variable", tstr); + } } /* ---------------------------------------------------------------------- */ @@ -467,7 +479,11 @@ void FixWallGran::post_force(int /*vflag*/) model->radj = 0.0; model->vj = vwall; model->omegaj = w0; - if (heat_flag) model->Tj = Twall; + if (heat_flag) { + if (tstr) + Twall = input->variable->compute_equal(tvar); + model->Tj = Twall; + } for (int i = 0; i < nlocal; i++) { if (! mask[i] & groupbit) continue; @@ -509,7 +525,6 @@ void FixWallGran::post_force(int /*vflag*/) } // Reset model and copy initial geometric data - model->reset_contact(); model->dx[0] = dx; model->dx[1] = dy; model->dx[2] = dz; diff --git a/src/GRANULAR/fix_wall_gran.h b/src/GRANULAR/fix_wall_gran.h index b61ac24d54..fd4a0a0ec8 100644 --- a/src/GRANULAR/fix_wall_gran.h +++ b/src/GRANULAR/fix_wall_gran.h @@ -67,6 +67,9 @@ class FixWallGran : public Fix { int heat_flag; int limit_damping; + int tvar; + char *tstr; + // shear history for single contact per particle double **history_one; diff --git a/src/GRANULAR/fix_wall_gran_region.cpp b/src/GRANULAR/fix_wall_gran_region.cpp index 7a774132f0..155b872e10 100644 --- a/src/GRANULAR/fix_wall_gran_region.cpp +++ b/src/GRANULAR/fix_wall_gran_region.cpp @@ -12,7 +12,7 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing authors: Dan Bolintineanu (SNL) + Contributing authors: Dan Bolintineanu (SNL), Joel Clemmer (SNL) ------------------------------------------------------------------------- */ #include "fix_wall_gran_region.h" @@ -50,7 +50,7 @@ FixWallGranRegion::FixWallGranRegion(LAMMPS *lmp, int narg, char **arg) : tmax = region->tmax; c2r = new int[tmax]; - model->wall_type = RDUPLICATE; + model->contact_type = WALLREGION; // re-allocate atom-based arrays with nshear // do not register with Atom class, since parent class did that @@ -226,7 +226,6 @@ void FixWallGranRegion::post_force(int /*vflag*/) // Reset model and copy initial geometric data // A bit unncessary since region->contact[ic] stores r - model->reset_contact(); model->dx[0] = region->contact[ic].delx; model->dx[1] = region->contact[ic].dely; model->dx[2] = region->contact[ic].delz; diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 52e4a84ded..b376bb64c1 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -14,8 +14,8 @@ /* ---------------------------------------------------------------------- Contributing authors: - Dan Bolintineanu (SNL), Ishan Srivastava (SNL), Jeremy Lechman(SNL) - Leo Silbert (SNL), Gary Grest (SNL) + Dan Bolintineanu (SNL), Joel Clemmer (SNL), Ishan Srivastava (SNL), + Jeremy Lechman(SNL), Leo Silbert (SNL), Gary Grest (SNL) ----------------------------------------------------------------------- */ #include "pair_granular.h" @@ -196,7 +196,6 @@ void PairGranular::compute(int eflag, int vflag) jtype = type[j]; // Reset model and copy initial geometric data - models[itype][jtype]->reset_contact(); models[itype][jtype]->xi = x[i]; models[itype][jtype]->xj = x[j]; models[itype][jtype]->radi = radius[i]; @@ -346,31 +345,30 @@ void PairGranular::coeff(int narg, char **arg) // Construct new model within vector vec_models.emplace_back(lmp); - //Parse mandatory normal and tangential specifications + //Parse mandatory specification int iarg = 2; vec_models.back().init_model(std::string(arg[iarg++]), NORMAL); iarg = vec_models.back().normal_model->parse_coeffs(arg, iarg, narg); - if (strcmp(arg[iarg++], "tangential") == 0) { - if (iarg >= narg) - error->all(FLERR,"Illegal pair_coeff command, must specify " - "tangential model after tangential keyword"); - vec_models.back().init_model(std::string(arg[iarg++]), TANGENTIAL); - iarg = vec_models.back().tangential_model->parse_coeffs(arg, iarg, narg); - } else{ - error->all(FLERR, "Illegal pair_coeff command, 'tangential' keyword expected"); - } - //Parse optional arguments while (iarg < narg) { - if (strcmp(arg[iarg], "damping") == 0) { + + if (strcmp(arg[iarg++], "tangential") == 0) { + if (iarg >= narg) + error->all(FLERR,"Illegal pair_coeff command, must specify " + "tangential model after tangential keyword"); + vec_models.back().init_model(std::string(arg[iarg++]), TANGENTIAL); + iarg = vec_models.back().tangential_model->parse_coeffs(arg, iarg, narg); + + } else if (strcmp(arg[iarg], "damping") == 0) { iarg++; if (iarg >= narg) error->all(FLERR, "Illegal pair_coeff command, must specify " "damping model after damping keyword"); vec_models.back().init_model(std::string(arg[iarg++]), DAMPING); iarg = vec_models.back().damping_model->parse_coeffs(arg, iarg, narg); + } else if (strcmp(arg[iarg], "rolling") == 0) { iarg++; if (iarg >= narg) @@ -386,6 +384,7 @@ void PairGranular::coeff(int narg, char **arg) "twisting model after twisting keyword"); vec_models.back().init_model(std::string(arg[iarg++]), TWISTING); iarg = vec_models.back().twisting_model->parse_coeffs(arg, iarg, narg); + } else if (strcmp(arg[iarg], "heat") == 0) { iarg++; if (iarg >= narg) @@ -394,6 +393,7 @@ void PairGranular::coeff(int narg, char **arg) vec_models.back().init_model(std::string(arg[iarg++]), HEAT); iarg = vec_models.back().heat_model->parse_coeffs(arg, iarg, narg); heat_flag = 1; + } else if (strcmp(arg[iarg], "cutoff") == 0) { if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters for cutoff keyword"); @@ -576,6 +576,8 @@ double PairGranular::init_one(int i, int j) for (int k = 0; k < NSUBMODELS; k++) vec_models.back().sub_models[k]->history_index = models[i][i]->sub_models[k]->history_index; + + cutoff_type[i][j] = cutoff_type[j][i] = MAX(cutoff_type[i][i], cutoff_type[j][j]); } // Check if heat model is defined for all type combinations @@ -732,7 +734,6 @@ double PairGranular::single(int i, int j, int itype, int jtype, double *radius = atom->radius; // Reset model and copy initial geometric data - models[itype][jtype]->reset_contact(); models[itype][jtype]->xi = x[i]; models[itype][jtype]->xj = x[j]; models[itype][jtype]->radi = radius[i];