tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Adding variable temperature to fix gran/wall, misc fixes/updates

Browse Source
This commit is contained in:
jtclemm
2022-09-14 21:40:00 -06:00
parent 06953bd67a
commit 038f4a5210
11 changed files with 98 additions and 81 deletions
Download Patch File Download Diff File Expand all files Collapse all files

62
src/GRANULAR/contact.cpp
View File

@ -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) {

12
src/GRANULAR/contact.h
View File

@ -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

2
src/GRANULAR/contact_normal_models.cpp
View File

@ -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");
}

12
src/GRANULAR/contact_sub_models.cpp
View File

@ -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"

2
src/GRANULAR/contact_tangential_models.cpp
View File

@ -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);
}

16
src/GRANULAR/contact_twisting_models.cpp
View File

@ -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;

3
src/GRANULAR/contact_twisting_models.h
View File

@ -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;
};
/* ---------------------------------------------------------------------- */

29
src/GRANULAR/fix_wall_gran.cpp
View File

@ -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;

3
src/GRANULAR/fix_wall_gran.h
View File

@ -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;

5
src/GRANULAR/fix_wall_gran_region.cpp
View File

@ -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;

33
src/GRANULAR/pair_granular.cpp
View File

@ -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];