From dc119da5855c582c866bd0804350c03c23507749 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 28 Jun 2013 16:14:30 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@10148 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/GRANULAR/pair_gran_hertz_history.cpp | 32 ++++++++------ src/GRANULAR/pair_gran_hooke.cpp | 33 +++++++++------ src/GRANULAR/pair_gran_hooke_history.cpp | 54 ++++++++++++++++++------ src/GRANULAR/pair_gran_hooke_history.h | 10 +++-- src/RIGID/fix_rigid_small.cpp | 47 ++++++++++++++++++--- src/RIGID/fix_rigid_small.h | 6 +++ 6 files changed, 135 insertions(+), 47 deletions(-) diff --git a/src/GRANULAR/pair_gran_hertz_history.cpp b/src/GRANULAR/pair_gran_hertz_history.cpp index 66f00be97f..a2248ec241 100644 --- a/src/GRANULAR/pair_gran_hertz_history.cpp +++ b/src/GRANULAR/pair_gran_hertz_history.cpp @@ -27,6 +27,7 @@ #include "neighbor.h" #include "neigh_list.h" #include "comm.h" +#include "memory.h" #include "error.h" using namespace LAMMPS_NS; @@ -60,12 +61,23 @@ void PairGranHertzHistory::compute(int eflag, int vflag) int shearupdate = 1; if (update->setupflag) shearupdate = 0; - // update rigid body ptrs and values for ghost atoms if using FixRigid masses + // update rigid body info for owned & ghost atoms if using FixRigid masses + // body[i] = which body atom I is in, -1 if none + // mass_body = mass of each rigid body if (fix_rigid && neighbor->ago == 0) { int tmp; - body = (int *) fix_rigid->extract("body",tmp); - mass_rigid = (double *) fix_rigid->extract("masstotal",tmp); + int *body = (int *) fix_rigid->extract("body",tmp); + double *mass_body = (double *) fix_rigid->extract("masstotal",tmp); + if (atom->nmax > nmax) { + memory->destroy(mass_rigid); + nmax = atom->nmax; + memory->create(mass_rigid,nmax,"pair:mass_rigid"); + } + int nlocal = atom->nlocal; + for (i = 0; i < nlocal; i++) + if (body[i] >= 0) mass_rigid[i] = mass_body[body[i]]; + else mass_rigid[i] = 0.0; comm->forward_comm_pair(this); } @@ -164,8 +176,8 @@ void PairGranHertzHistory::compute(int eflag, int vflag) mj = mass[type[j]]; } if (fix_rigid) { - if (body[i] >= 0) mi = mass_rigid[body[i]]; - if (body[j] >= 0) mj = mass_rigid[body[j]]; + if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; + if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; } meff = mi*mj / (mi+mj); @@ -374,13 +386,9 @@ double PairGranHertzHistory::single(int i, int j, int itype, int jtype, mj = mass[type[j]]; } if (fix_rigid) { - // NOTE: need to make sure ghost atoms have updated body? - // depends on where single() is called from - int tmp; - body = (int *) fix_rigid->extract("body",tmp); - mass_rigid = (double *) fix_rigid->extract("masstotal",tmp); - if (body[i] >= 0) mi = mass_rigid[body[i]]; - if (body[j] >= 0) mj = mass_rigid[body[j]]; + // NOTE: insure mass_rigid is current for owned+ghost atoms? + if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; + if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; } meff = mi*mj / (mi+mj); diff --git a/src/GRANULAR/pair_gran_hooke.cpp b/src/GRANULAR/pair_gran_hooke.cpp index 9272909870..c6b1a83231 100644 --- a/src/GRANULAR/pair_gran_hooke.cpp +++ b/src/GRANULAR/pair_gran_hooke.cpp @@ -25,6 +25,7 @@ #include "neighbor.h" #include "neigh_list.h" #include "comm.h" +#include "memory.h" using namespace LAMMPS_NS; @@ -53,12 +54,23 @@ void PairGranHooke::compute(int eflag, int vflag) if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = 0; - // update rigid body ptrs and values for ghost atoms if using FixRigid masses + // update rigid body info for owned & ghost atoms if using FixRigid masses + // body[i] = which body atom I is in, -1 if none + // mass_body = mass of each rigid body if (fix_rigid && neighbor->ago == 0) { int tmp; - body = (int *) fix_rigid->extract("body",tmp); - mass_rigid = (double *) fix_rigid->extract("masstotal",tmp); + int *body = (int *) fix_rigid->extract("body",tmp); + double *mass_body = (double *) fix_rigid->extract("masstotal",tmp); + if (atom->nmax > nmax) { + memory->destroy(mass_rigid); + nmax = atom->nmax; + memory->create(mass_rigid,nmax,"pair:mass_rigid"); + } + int nlocal = atom->nlocal; + for (i = 0; i < nlocal; i++) + if (body[i] >= 0) mass_rigid[i] = mass_body[body[i]]; + else mass_rigid[i] = 0.0; comm->forward_comm_pair(this); } @@ -144,8 +156,8 @@ void PairGranHooke::compute(int eflag, int vflag) mj = mass[type[j]]; } if (fix_rigid) { - if (body[i] >= 0) mi = mass_rigid[body[i]]; - if (body[j] >= 0) mj = mass_rigid[body[j]]; + if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; + if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; } meff = mi*mj / (mi+mj); @@ -290,14 +302,11 @@ double PairGranHooke::single(int i, int j, int itype, int jtype, double rsq, mi = mass[type[i]]; mj = mass[type[j]]; } + if (fix_rigid) { - // NOTE: need to make sure ghost atoms have updated body? - // depends on where single() is called from - int tmp; - body = (int *) fix_rigid->extract("body",tmp); - mass_rigid = (double *) fix_rigid->extract("masstotal",tmp); - if (body[i] >= 0) mi = mass_rigid[body[i]]; - if (body[j] >= 0) mj = mass_rigid[body[j]]; + // NOTE: insure mass_rigid is current for owned+ghost atoms? + if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; + if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; } meff = mi*mj / (mi+mj); diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp index 115d19f602..ec9b2e559d 100644 --- a/src/GRANULAR/pair_gran_hooke_history.cpp +++ b/src/GRANULAR/pair_gran_hooke_history.cpp @@ -53,6 +53,13 @@ PairGranHookeHistory::PairGranHookeHistory(LAMMPS *lmp) : Pair(lmp) computeflag = 0; neighprev = 0; + + nmax = 0; + mass_rigid = NULL; + + // set comm size needed by this Pair if used with fix rigid + + comm_forward = 1; } /* ---------------------------------------------------------------------- */ @@ -72,6 +79,8 @@ PairGranHookeHistory::~PairGranHookeHistory() delete [] maxrad_dynamic; delete [] maxrad_frozen; } + + memory->destroy(mass_rigid); } /* ---------------------------------------------------------------------- */ @@ -98,12 +107,23 @@ void PairGranHookeHistory::compute(int eflag, int vflag) int shearupdate = 1; if (update->setupflag) shearupdate = 0; - // update rigid body ptrs and values for ghost atoms if using FixRigid masses + // update rigid body info for owned & ghost atoms if using FixRigid masses + // body[i] = which body atom I is in, -1 if none + // mass_body = mass of each rigid body if (fix_rigid && neighbor->ago == 0) { int tmp; - body = (int *) fix_rigid->extract("body",tmp); - mass_rigid = (double *) fix_rigid->extract("masstotal",tmp); + int *body = (int *) fix_rigid->extract("body",tmp); + double *mass_body = (double *) fix_rigid->extract("masstotal",tmp); + if (atom->nmax > nmax) { + memory->destroy(mass_rigid); + nmax = atom->nmax; + memory->create(mass_rigid,nmax,"pair:mass_rigid"); + } + int nlocal = atom->nlocal; + for (i = 0; i < nlocal; i++) + if (body[i] >= 0) mass_rigid[i] = mass_body[body[i]]; + else mass_rigid[i] = 0.0; comm->forward_comm_pair(this); } @@ -202,8 +222,8 @@ void PairGranHookeHistory::compute(int eflag, int vflag) mj = mass[type[j]]; } if (fix_rigid) { - if (body[i] >= 0) mi = mass_rigid[body[i]]; - if (body[j] >= 0) mj = mass_rigid[body[j]]; + if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; + if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; } meff = mi*mj / (mi+mj); @@ -653,13 +673,9 @@ double PairGranHookeHistory::single(int i, int j, int itype, int jtype, mj = mass[type[j]]; } if (fix_rigid) { - // NOTE: need to make sure ghost atoms have updated body? - // depends on where single() is called from - int tmp; - body = (int *) fix_rigid->extract("body",tmp); - mass_rigid = (double *) fix_rigid->extract("masstotal",tmp); - if (body[i] >= 0) mi = mass_rigid[body[i]]; - if (body[j] >= 0) mj = mass_rigid[body[j]]; + // NOTE: insure mass_rigid is current for owned+ghost atoms? + if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; + if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; } meff = mi*mj / (mi+mj); @@ -745,7 +761,7 @@ int PairGranHookeHistory::pack_comm(int n, int *list, m = 0; for (i = 0; i < n; i++) { j = list[i]; - buf[m++] = body[j]; + buf[m++] = mass_rigid[j]; } return 1; } @@ -759,7 +775,7 @@ void PairGranHookeHistory::unpack_comm(int n, int first, double *buf) m = 0; last = first + n; for (i = first; i < last; i++) - body[i] = static_cast (buf[m++]); + mass_rigid[i] = buf[m++]; } /* ---------------------------------------------------------------------- */ @@ -770,3 +786,13 @@ void *PairGranHookeHistory::extract(const char *str, int &dim) if (strcmp(str,"computeflag") == 0) return (void *) &computeflag; return NULL; } + +/* ---------------------------------------------------------------------- + memory usage of local atom-based arrays +------------------------------------------------------------------------- */ + +double PairGranHookeHistory::memory_usage() +{ + double bytes = nmax * sizeof(double); + return bytes; +} diff --git a/src/GRANULAR/pair_gran_hooke_history.h b/src/GRANULAR/pair_gran_hooke_history.h index 06d720ea04..ce69e382c6 100644 --- a/src/GRANULAR/pair_gran_hooke_history.h +++ b/src/GRANULAR/pair_gran_hooke_history.h @@ -45,6 +45,7 @@ class PairGranHookeHistory : public Pair { int pack_comm(int, int *, double *, int, int *); void unpack_comm(int, int, double *); void *extract(const char *, int &); + double memory_usage(); protected: double kn,kt,gamman,gammat,xmu; @@ -59,9 +60,12 @@ class PairGranHookeHistory : public Pair { double *maxrad_dynamic,*maxrad_frozen; class FixShearHistory *fix_history; - class Fix *fix_rigid; - int *body; - double *mass_rigid; + + // storage of rigid body masses for use in granular interactions + + class Fix *fix_rigid; // ptr to rigid body fix, NULL if none + double *mass_rigid; // rigid mass for owned+ghost atoms + int nmax; // allocated size of mass_rigid void allocate(); }; diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index dba76f60da..31190dcd9d 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/src/RIGID/fix_rigid_small.cpp @@ -151,6 +151,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : bodyown[i] = nlocal_body++; } else bodyown[i] = -1; + // bodysize = sizeof(Body) in doubles bodysize = sizeof(Body)/sizeof(double); @@ -217,6 +218,11 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : random = NULL; if (langflag) random = new RanMars(lmp,seed + comm->me); + // mass vector for granular pair styles + + mass_body = NULL; + nmax_mass = 0; + // firstflag = 1 triggers one-time initialization of rigid body attributes firstflag = 1; @@ -244,6 +250,7 @@ FixRigidSmall::~FixRigidSmall() delete random; memory->destroy(langextra); + memory->destroy(mass_body); } /* ---------------------------------------------------------------------- */ @@ -1289,10 +1296,7 @@ void FixRigidSmall::create_bodies() n = 0; for (i = 0; i < nlocal; i++) { if (!(mask[i] & groupbit)) continue; - if (hash->find(molecule[i]) == hash->end()) { - hash->insert(std::pair (molecule[i],n)); - n++; - } + if (hash->find(molecule[i]) == hash->end()) (*hash)[molecule[i]] = n++; } // bbox = bounding box of each rigid body my atoms are part of @@ -1552,7 +1556,7 @@ void FixRigidSmall::setup_bodies() atom->tri_flag || atom->mu_flag) { int flag = 0; for (i = 0; i < nlocal; i++) { - if (atom2body[i] < 0) continue; + if (bodytag[i] == 0) continue; if (radius && radius[i] > 0.0) flag = 1; if (ellipsoid && ellipsoid[i] >= 0) flag = 1; if (line && line[i] >= 0) flag = 1; @@ -1577,7 +1581,7 @@ void FixRigidSmall::setup_bodies() for (i = 0; i < nlocal; i++) { eflags[i] = 0; - if (atom2body[i] < 0) continue; + if (bodytag[i] == 0) continue; // set to POINT or SPHERE or ELLIPSOID or LINE @@ -2518,6 +2522,37 @@ void FixRigidSmall::reset_dt() dtq = 0.5 * update->dt; } +/* ---------------------------------------------------------------------- */ + +void *FixRigidSmall::extract(const char *str, int &dim) +{ + if (strcmp(str,"body") == 0) { + dim = 1; + return atom2body; + } + + // return vector of rigid body masses, for owned+ghost bodies + // used by granular pair styles, indexed by atom2body + + if (strcmp(str,"masstotal") == 0) { + dim = 1; + + if (nmax_mass < nmax_body) { + memory->destroy(mass_body); + nmax_mass = nmax_body; + memory->create(mass_body,nmax_mass,"rigid:mass_body"); + } + + int n = nlocal_body + nghost_body; + for (int i = 0; i < n; i++) + mass_body[i] = body[i].mass; + + return mass_body; + } + + return NULL; +} + /* ---------------------------------------------------------------------- return temperature of collection of rigid bodies non-active DOF are removed by fflag/tflag and in tfactor diff --git a/src/RIGID/fix_rigid_small.h b/src/RIGID/fix_rigid_small.h index 8eca0e8d5e..db5d0adb2f 100644 --- a/src/RIGID/fix_rigid_small.h +++ b/src/RIGID/fix_rigid_small.h @@ -59,6 +59,7 @@ class FixRigidSmall : public Fix { int dof(int); void deform(int); void reset_dt(); + void *extract(const char*,int &); double compute_scalar(); double memory_usage(); @@ -128,6 +129,11 @@ class FixRigidSmall : public Fix { int **counts; // counts of atom types in bodies double **itensor; // 6 space-frame components of inertia tensor + // mass per body, accessed by granular pair styles + + double *mass_body; + int nmax_mass; + // Langevin thermostatting int langflag; // 0/1 = no/yes Langevin thermostat