diff --git a/src/SPIN/atom_vec_spin.cpp b/src/SPIN/atom_vec_spin.cpp index a8ff15492a..638a3b8021 100644 --- a/src/SPIN/atom_vec_spin.cpp +++ b/src/SPIN/atom_vec_spin.cpp @@ -62,14 +62,16 @@ AtomVecSpin::AtomVecSpin(LAMMPS *lmp) : AtomVec(lmp) } /* ---------------------------------------------------------------------- - clear all forces (mechanical and magnetic) + clear extra forces starting at atom N + nbytes = # of bytes to clear for a per-atom vector + include f b/c this is invoked from within SPIN pair styles ------------------------------------------------------------------------- */ -void AtomVecSpin::force_clear(int /*n*/, size_t nbytes) +void AtomVecSpin::force_clear(int n, size_t nbytes) { - memset(&atom->f[0][0],0,3*nbytes); - memset(&atom->fm[0][0],0,3*nbytes); - memset(&atom->fm_long[0][0],0,3*nbytes); + memset(&atom->f[n][0],0,3*nbytes); + memset(&atom->fm[n][0],0,3*nbytes); + memset(&atom->fm_long[n][0],0,3*nbytes); } /* ---------------------------------------------------------------------- diff --git a/src/USER-EFF/atom_vec_electron.cpp b/src/USER-EFF/atom_vec_electron.cpp index 0bad0f115f..3688f7f582 100644 --- a/src/USER-EFF/atom_vec_electron.cpp +++ b/src/USER-EFF/atom_vec_electron.cpp @@ -18,14 +18,8 @@ #include "atom_vec_electron.h" #include #include "atom.h" -#include "comm.h" -#include "domain.h" -#include "modify.h" -#include "fix.h" #include "citeme.h" -#include "memory.h" #include "error.h" -#include "utils.h" using namespace LAMMPS_NS; @@ -46,6 +40,8 @@ AtomVecElectron::AtomVecElectron(LAMMPS *lmp) : AtomVec(lmp) { if (lmp->citeme) lmp->citeme->add(cite_user_eff_package); + mass_type = 1; + molecular = 0; forceclearflag = 1; atom->ecp_flag = 0; @@ -75,11 +71,14 @@ AtomVecElectron::AtomVecElectron(LAMMPS *lmp) : AtomVec(lmp) setup_fields(); } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + clear extra forces starting at atom N + nbytes = # of bytes to clear for a per-atom vector +------------------------------------------------------------------------- */ -void AtomVecElectron::force_clear(int /*n*/, size_t nbytes) +void AtomVecElectron::force_clear(int n, size_t nbytes) { - memset(&atom->erforce[0],0,nbytes); + memset(&atom->erforce[n],0,nbytes); } /* ---------------------------------------------------------------------- diff --git a/src/USER-MESO/atom_vec_edpd.cpp b/src/USER-MESO/atom_vec_edpd.cpp index b4602c760c..7a031c256b 100644 --- a/src/USER-MESO/atom_vec_edpd.cpp +++ b/src/USER-MESO/atom_vec_edpd.cpp @@ -29,9 +29,6 @@ using namespace LAMMPS_NS; AtomVecEDPD::AtomVecEDPD(LAMMPS *lmp) : AtomVec(lmp) { - if (strcmp(update->unit_style,"lj") != 0) - error->all(FLERR,"Atom style edpd requires lj units"); - molecular = 0; mass_type = 1; forceclearflag = 1; @@ -62,9 +59,22 @@ AtomVecEDPD::AtomVecEDPD(LAMMPS *lmp) : AtomVec(lmp) /* ---------------------------------------------------------------------- */ -void AtomVecEDPD::force_clear(int /*n*/, size_t nbytes) +void AtomVecEDPD::init() { - memset(&atom->edpd_flux[0],0,nbytes); + AtomVec::init(); + + if (strcmp(update->unit_style,"lj") != 0) + error->all(FLERR,"Atom style edpd requires lj units"); +} + +/* ---------------------------------------------------------------------- + clear extra forces starting at atom N + nbytes = # of bytes to clear for a per-atom vector +------------------------------------------------------------------------- */ + +void AtomVecEDPD::force_clear(int n, size_t nbytes) +{ + memset(&atom->edpd_flux[n],0,nbytes); } /* ---------------------------------------------------------------------- diff --git a/src/USER-MESO/atom_vec_edpd.h b/src/USER-MESO/atom_vec_edpd.h index 7d41b51665..bb667ae792 100644 --- a/src/USER-MESO/atom_vec_edpd.h +++ b/src/USER-MESO/atom_vec_edpd.h @@ -27,6 +27,7 @@ namespace LAMMPS_NS { class AtomVecEDPD : public AtomVec { public: AtomVecEDPD(class LAMMPS *); + void init(); void force_clear(int, size_t); void create_atom_post(int); void data_atom_post(int); diff --git a/src/USER-MESO/atom_vec_mdpd.cpp b/src/USER-MESO/atom_vec_mdpd.cpp index 4c9db36645..eefda1ff6a 100644 --- a/src/USER-MESO/atom_vec_mdpd.cpp +++ b/src/USER-MESO/atom_vec_mdpd.cpp @@ -14,14 +14,8 @@ #include "atom_vec_mdpd.h" #include #include "atom.h" -#include "comm.h" -#include "domain.h" -#include "modify.h" -#include "fix.h" #include "update.h" -#include "memory.h" #include "error.h" -#include "utils.h" using namespace LAMMPS_NS; @@ -29,865 +23,63 @@ using namespace LAMMPS_NS; AtomVecMDPD::AtomVecMDPD(LAMMPS *lmp) : AtomVec(lmp) { - if(strcmp(update->unit_style,"lj") != 0) - error->all(FLERR,"Atom style mdpd requires lj units"); - molecular = 0; mass_type = 1; forceclearflag = 1; - comm_x_only = comm_f_only = 0; - comm->ghost_velocity = 1; - - size_forward = 3 + 4; // 3 + rho + vest[3], that means we may only communicate 4 in hybrid - size_reverse = 3 + 1; // 3 + drho - size_border = 6 + 4; // 6 + rho + vest[3] - size_velocity = 3; - size_data_atom = 6; - size_data_vel = 4; - xcol_data = 4; - atom->rho_flag = 1; atom->vest_flag = 1; + + // strings with peratom variables to include in each AtomVec method + // strings cannot contain fields in corresponding AtomVec default strings + // order of fields in a string does not matter + // except: fields_data_atom & fields_data_vel must match data file + + fields_grow = (char *) "rho drho vest"; + fields_copy = (char *) "rho drho vest"; + fields_comm = (char *) "rho vest"; + fields_comm_vel = (char *) "rho vest"; + fields_reverse = (char *) "drho"; + fields_border = (char *) "rho vest"; + fields_border_vel = (char *) "rho vest"; + fields_exchange = (char *) "rho vest"; + fields_restart = (char * ) "rho vest"; + fields_create = (char *) "rho drho vest"; + fields_data_atom = (char *) "id type rho x"; + fields_data_vel = (char *) "id v"; + + setup_fields(); } -/* ---------------------------------------------------------------------- - grow atom arrays - n = 0 grows arrays by a chunk - n > 0 allocates arrays to size n - ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- */ -void AtomVecMDPD::grow(int n) +void AtomVecMDPD::init() { - if (n == 0) grow_nmax(); - else nmax = n; - atom->nmax = nmax; - if (nmax < 0 || nmax > MAXSMALLINT) - error->one(FLERR,"Per-processor system is too big"); - - tag = memory->grow(atom->tag, nmax, "atom:tag"); - type = memory->grow(atom->type, nmax, "atom:type"); - mask = memory->grow(atom->mask, nmax, "atom:mask"); - image = memory->grow(atom->image, nmax, "atom:image"); - x = memory->grow(atom->x, nmax, 3, "atom:x"); - v = memory->grow(atom->v, nmax, 3, "atom:v"); - f = memory->grow(atom->f, nmax*comm->nthreads, 3, "atom:f"); - - rho = memory->grow(atom->rho, nmax, "atom:rho"); - drho = memory->grow(atom->drho, nmax*comm->nthreads, "atom:drho"); - vest = memory->grow(atom->vest, nmax, 3, "atom:vest"); - - if (atom->nextra_grow) - for (int iextra = 0; iextra < atom->nextra_grow; iextra++) - modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax); + if (strcmp(update->unit_style,"lj") != 0) + error->all(FLERR,"Atom style mdpd requires lj units"); } /* ---------------------------------------------------------------------- - reset local array ptrs - ------------------------------------------------------------------------- */ - -void AtomVecMDPD::grow_reset() { - tag = atom->tag; - type = atom->type; - mask = atom->mask; - image = atom->image; - x = atom->x; - v = atom->v; - f = atom->f; - rho = atom->rho; - drho = atom->drho; - vest = atom->vest; -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecMDPD::copy(int i, int j, int delflag) { - //printf("in AtomVecMDPD::copy\n"); - tag[j] = tag[i]; - type[j] = type[i]; - mask[j] = mask[i]; - image[j] = image[i]; - x[j][0] = x[i][0]; - x[j][1] = x[i][1]; - x[j][2] = x[i][2]; - v[j][0] = v[i][0]; - v[j][1] = v[i][1]; - v[j][2] = v[i][2]; - - rho[j] = rho[i]; - drho[j] = drho[i]; - vest[j][0] = vest[i][0]; - vest[j][1] = vest[i][1]; - vest[j][2] = vest[i][2]; - - if (atom->nextra_grow) - for (int iextra = 0; iextra < atom->nextra_grow; iextra++) - modify->fix[atom->extra_grow[iextra]]->copy_arrays(i, j,delflag); -} - -/* ---------------------------------------------------------------------- */ + clear extra forces starting at atom N + nbytes = # of bytes to clear for a per-atom vector +------------------------------------------------------------------------- */ void AtomVecMDPD::force_clear(int n, size_t nbytes) { - memset(&drho[n],0,nbytes); -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecMDPD::pack_comm_hybrid(int n, int *list, double *buf) { - //printf("in AtomVecMDPD::pack_comm_hybrid\n"); - int i, j, m; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = rho[j]; - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecMDPD::unpack_comm_hybrid(int n, int first, double *buf) { - //printf("in AtomVecMDPD::unpack_comm_hybrid\n"); - int i, m, last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - rho[i] = buf[m++]; - vest[i][0] = buf[m++]; - vest[i][1] = buf[m++]; - vest[i][2] = buf[m++]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecMDPD::pack_border_hybrid(int n, int *list, double *buf) { - //printf("in AtomVecMDPD::pack_border_hybrid\n"); - int i, j, m; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = rho[j]; - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecMDPD::unpack_border_hybrid(int n, int first, double *buf) { - //printf("in AtomVecMDPD::unpack_border_hybrid\n"); - int i, m, last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - rho[i] = buf[m++]; - vest[i][0] = buf[m++]; - vest[i][1] = buf[m++]; - vest[i][2] = buf[m++]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecMDPD::pack_reverse_hybrid(int n, int first, double *buf) { - //printf("in AtomVecMDPD::pack_reverse_hybrid\n"); - int i, m, last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - buf[m++] = drho[i]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecMDPD::unpack_reverse_hybrid(int n, int *list, double *buf) { - //printf("in AtomVecMDPD::unpack_reverse_hybrid\n"); - int i, j, m; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - drho[j] += buf[m++]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecMDPD::pack_comm(int n, int *list, double *buf, int pbc_flag, - int *pbc) { - //printf("in AtomVecMDPD::pack_comm\n"); - int i, j, m; - double dx, dy, dz; - - m = 0; - if (pbc_flag == 0) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0]; - buf[m++] = x[j][1]; - buf[m++] = x[j][2]; - buf[m++] = rho[j]; - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } else { - if (domain->triclinic == 0) { - dx = pbc[0] * domain->xprd; - dy = pbc[1] * domain->yprd; - dz = pbc[2] * domain->zprd; - } else { - dx = pbc[0] * domain->xprd + pbc[5] * domain->xy + pbc[4] * domain->xz; - dy = pbc[1] * domain->yprd + pbc[3] * domain->yz; - dz = pbc[2] * domain->zprd; - } - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0] + dx; - buf[m++] = x[j][1] + dy; - buf[m++] = x[j][2] + dz; - buf[m++] = rho[j]; - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecMDPD::pack_comm_vel(int n, int *list, double *buf, int pbc_flag, - int *pbc) { - //printf("in AtomVecMDPD::pack_comm_vel\n"); - int i, j, m; - double dx, dy, dz; - - m = 0; - if (pbc_flag == 0) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0]; - buf[m++] = x[j][1]; - buf[m++] = x[j][2]; - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; - buf[m++] = rho[j]; - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } else { - if (domain->triclinic == 0) { - dx = pbc[0] * domain->xprd; - dy = pbc[1] * domain->yprd; - dz = pbc[2] * domain->zprd; - } else { - dx = pbc[0] * domain->xprd + pbc[5] * domain->xy + pbc[4] * domain->xz; - dy = pbc[1] * domain->yprd + pbc[3] * domain->yz; - dz = pbc[2] * domain->zprd; - } - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0] + dx; - buf[m++] = x[j][1] + dy; - buf[m++] = x[j][2] + dz; - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; - buf[m++] = rho[j]; - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecMDPD::unpack_comm(int n, int first, double *buf) { - //printf("in AtomVecMDPD::unpack_comm\n"); - int i, m, last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - x[i][0] = buf[m++]; - x[i][1] = buf[m++]; - x[i][2] = buf[m++]; - rho[i] = buf[m++]; - vest[i][0] = buf[m++]; - vest[i][1] = buf[m++]; - vest[i][2] = buf[m++]; - } -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecMDPD::unpack_comm_vel(int n, int first, double *buf) { - //printf("in AtomVecMDPD::unpack_comm_vel\n"); - int i, m, last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - x[i][0] = buf[m++]; - x[i][1] = buf[m++]; - x[i][2] = buf[m++]; - v[i][0] = buf[m++]; - v[i][1] = buf[m++]; - v[i][2] = buf[m++]; - rho[i] = buf[m++]; - vest[i][0] = buf[m++]; - vest[i][1] = buf[m++]; - vest[i][2] = buf[m++]; - } -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecMDPD::pack_reverse(int n, int first, double *buf) { - //printf("in AtomVecMDPD::pack_reverse\n"); - int i, m, last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - buf[m++] = f[i][0]; - buf[m++] = f[i][1]; - buf[m++] = f[i][2]; - buf[m++] = drho[i]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecMDPD::unpack_reverse(int n, int *list, double *buf) { - //printf("in AtomVecMDPD::unpack_reverse\n"); - int i, j, m; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - f[j][0] += buf[m++]; - f[j][1] += buf[m++]; - f[j][2] += buf[m++]; - drho[j] += buf[m++]; - } -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecMDPD::pack_border(int n, int *list, double *buf, int pbc_flag, - int *pbc) { - //printf("in AtomVecMDPD::pack_border\n"); - int i, j, m; - double dx, dy, dz; - - m = 0; - if (pbc_flag == 0) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0]; - buf[m++] = x[j][1]; - buf[m++] = x[j][2]; - buf[m++] = ubuf(tag[j]).d; - buf[m++] = ubuf(type[j]).d; - buf[m++] = ubuf(mask[j]).d; - buf[m++] = rho[j]; - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } else { - if (domain->triclinic == 0) { - dx = pbc[0] * domain->xprd; - dy = pbc[1] * domain->yprd; - dz = pbc[2] * domain->zprd; - } else { - dx = pbc[0]; - dy = pbc[1]; - dz = pbc[2]; - } - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0] + dx; - buf[m++] = x[j][1] + dy; - buf[m++] = x[j][2] + dz; - buf[m++] = ubuf(tag[j]).d; - buf[m++] = ubuf(type[j]).d; - buf[m++] = ubuf(mask[j]).d; - buf[m++] = rho[j]; - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } - - if (atom->nextra_border) - for (int iextra = 0; iextra < atom->nextra_border; iextra++) - m += modify->fix[atom->extra_border[iextra]]->pack_border(n,list,&buf[m]); - - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecMDPD::pack_border_vel(int n, int *list, double *buf, int pbc_flag, - int *pbc) -{ - int i,j,m; - double dx,dy,dz,dvx,dvy,dvz; - - m = 0; - if (pbc_flag == 0) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0]; - buf[m++] = x[j][1]; - buf[m++] = x[j][2]; - buf[m++] = ubuf(tag[j]).d; - buf[m++] = ubuf(type[j]).d; - buf[m++] = ubuf(mask[j]).d; - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; - buf[m++] = rho[j]; - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } else { - if (domain->triclinic == 0) { - dx = pbc[0] * domain->xprd; - dy = pbc[1] * domain->yprd; - dz = pbc[2] * domain->zprd; - } else { - dx = pbc[0]; - dy = pbc[1]; - dz = pbc[2]; - } - if (!deform_vremap) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0] + dx; - buf[m++] = x[j][1] + dy; - buf[m++] = x[j][2] + dz; - buf[m++] = ubuf(tag[j]).d; - buf[m++] = ubuf(type[j]).d; - buf[m++] = ubuf(mask[j]).d; - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; - buf[m++] = rho[j]; - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } else { - dvx = pbc[0] * h_rate[0] + pbc[5] * h_rate[5] + pbc[4] * h_rate[4]; - dvy = pbc[1] * h_rate[1] + pbc[3] * h_rate[3]; - dvz = pbc[2] * h_rate[2]; - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0] + dx; - buf[m++] = x[j][1] + dy; - buf[m++] = x[j][2] + dz; - buf[m++] = ubuf(tag[j]).d; - buf[m++] = ubuf(type[j]).d; - buf[m++] = ubuf(mask[j]).d; - if (mask[i] & deform_groupbit) { - buf[m++] = v[j][0] + dvx; - buf[m++] = v[j][1] + dvy; - buf[m++] = v[j][2] + dvz; - buf[m++] = rho[j]; - buf[m++] = vest[j][0] + dvx; - buf[m++] = vest[j][1] + dvy; - buf[m++] = vest[j][2] + dvz; - } else { - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; - buf[m++] = rho[j]; - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } - } - } - - if (atom->nextra_border) - for (int iextra = 0; iextra < atom->nextra_border; iextra++) - m += modify->fix[atom->extra_border[iextra]]->pack_border(n,list,&buf[m]); - - return m; -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecMDPD::unpack_border(int n, int first, double *buf) { - //printf("in AtomVecMDPD::unpack_border\n"); - int i, m, last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - if (i == nmax) - grow(0); - x[i][0] = buf[m++]; - x[i][1] = buf[m++]; - x[i][2] = buf[m++]; - tag[i] = (tagint) ubuf(buf[m++]).i; - type[i] = (int) ubuf(buf[m++]).i; - mask[i] = (int) ubuf(buf[m++]).i; - rho[i] = buf[m++]; - vest[i][0] = buf[m++]; - vest[i][1] = buf[m++]; - vest[i][2] = buf[m++]; - } - - if (atom->nextra_border) - for (int iextra = 0; iextra < atom->nextra_border; iextra++) - m += modify->fix[atom->extra_border[iextra]]-> - unpack_border(n,first,&buf[m]); -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecMDPD::unpack_border_vel(int n, int first, double *buf) { - //printf("in AtomVecMDPD::unpack_border_vel\n"); - int i, m, last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - if (i == nmax) - grow(0); - x[i][0] = buf[m++]; - x[i][1] = buf[m++]; - x[i][2] = buf[m++]; - tag[i] = (tagint) ubuf(buf[m++]).i; - type[i] = (int) ubuf(buf[m++]).i; - mask[i] = (int) ubuf(buf[m++]).i; - v[i][0] = buf[m++]; - v[i][1] = buf[m++]; - v[i][2] = buf[m++]; - rho[i] = buf[m++]; - vest[i][0] = buf[m++]; - vest[i][1] = buf[m++]; - vest[i][2] = buf[m++]; - } - - if (atom->nextra_border) - for (int iextra = 0; iextra < atom->nextra_border; iextra++) - m += modify->fix[atom->extra_border[iextra]]-> - unpack_border(n,first,&buf[m]); + memset(&atom->drho[n],0,nbytes); } /* ---------------------------------------------------------------------- - pack data for atom I for sending to another proc - xyz must be 1st 3 values, so comm::exchange() can test on them - ------------------------------------------------------------------------- */ - -int AtomVecMDPD::pack_exchange(int i, double *buf) { - //printf("in AtomVecMDPD::pack_exchange\n"); - int m = 1; - buf[m++] = x[i][0]; - buf[m++] = x[i][1]; - buf[m++] = x[i][2]; - buf[m++] = v[i][0]; - buf[m++] = v[i][1]; - buf[m++] = v[i][2]; - buf[m++] = ubuf(tag[i]).d; - buf[m++] = ubuf(type[i]).d; - buf[m++] = ubuf(mask[i]).d; - buf[m++] = ubuf(image[i]).d; - buf[m++] = rho[i]; - buf[m++] = vest[i][0]; - buf[m++] = vest[i][1]; - buf[m++] = vest[i][2]; - - if (atom->nextra_grow) - for (int iextra = 0; iextra < atom->nextra_grow; iextra++) - m += modify->fix[atom->extra_grow[iextra]]->pack_exchange(i, &buf[m]); - - buf[0] = m; - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecMDPD::unpack_exchange(double *buf) { - //printf("in AtomVecMDPD::unpack_exchange\n"); - int nlocal = atom->nlocal; - if (nlocal == nmax) - grow(0); - - int m = 1; - x[nlocal][0] = buf[m++]; - x[nlocal][1] = buf[m++]; - x[nlocal][2] = buf[m++]; - v[nlocal][0] = buf[m++]; - v[nlocal][1] = buf[m++]; - v[nlocal][2] = buf[m++]; - tag[nlocal] = (tagint) ubuf(buf[m++]).i; - type[nlocal] = (int) ubuf(buf[m++]).i; - mask[nlocal] = (int) ubuf(buf[m++]).i; - image[nlocal] = (imageint) ubuf(buf[m++]).i; - rho[nlocal] = buf[m++]; - vest[nlocal][0] = buf[m++]; - vest[nlocal][1] = buf[m++]; - vest[nlocal][2] = buf[m++]; - - if (atom->nextra_grow) - for (int iextra = 0; iextra < atom->nextra_grow; iextra++) - m += modify->fix[atom->extra_grow[iextra]]-> unpack_exchange(nlocal, - &buf[m]); - - atom->nlocal++; - return m; -} - -/* ---------------------------------------------------------------------- - size of restart data for all atoms owned by this proc - include extra data stored by fixes - ------------------------------------------------------------------------- */ - -int AtomVecMDPD::size_restart() { - int i; - - int nlocal = atom->nlocal; - int n = 15 * nlocal; // 11 + rho + vest[3] - - if (atom->nextra_restart) - for (int iextra = 0; iextra < atom->nextra_restart; iextra++) - for (i = 0; i < nlocal; i++) - n += modify->fix[atom->extra_restart[iextra]]->size_restart(i); - - return n; -} - -/* ---------------------------------------------------------------------- - pack atom I's data for restart file including extra quantities - xyz must be 1st 3 values, so that read_restart can test on them - molecular types may be negative, but write as positive - ------------------------------------------------------------------------- */ - -int AtomVecMDPD::pack_restart(int i, double *buf) { - int m = 1; - buf[m++] = x[i][0]; - buf[m++] = x[i][1]; - buf[m++] = x[i][2]; - buf[m++] = ubuf(tag[i]).d; - buf[m++] = ubuf(type[i]).d; - buf[m++] = ubuf(mask[i]).d; - buf[m++] = ubuf(image[i]).d; - buf[m++] = v[i][0]; - buf[m++] = v[i][1]; - buf[m++] = v[i][2]; - buf[m++] = rho[i]; - buf[m++] = vest[i][0]; - buf[m++] = vest[i][1]; - buf[m++] = vest[i][2]; - - if (atom->nextra_restart) - for (int iextra = 0; iextra < atom->nextra_restart; iextra++) - m += modify->fix[atom->extra_restart[iextra]]->pack_restart(i, &buf[m]); - - buf[0] = m; - return m; -} - -/* ---------------------------------------------------------------------- - unpack data for one atom from restart file including extra quantities - ------------------------------------------------------------------------- */ - -int AtomVecMDPD::unpack_restart(double *buf) { - int nlocal = atom->nlocal; - if (nlocal == nmax) { - grow(0); - if (atom->nextra_store) - memory->grow(atom->extra, nmax, atom->nextra_store, "atom:extra"); - } - - int m = 1; - x[nlocal][0] = buf[m++]; - x[nlocal][1] = buf[m++]; - x[nlocal][2] = buf[m++]; - tag[nlocal] = (tagint) ubuf(buf[m++]).i; - type[nlocal] = (int) ubuf(buf[m++]).i; - mask[nlocal] = (int) ubuf(buf[m++]).i; - image[nlocal] = (imageint) ubuf(buf[m++]).i; - v[nlocal][0] = buf[m++]; - v[nlocal][1] = buf[m++]; - v[nlocal][2] = buf[m++]; - rho[nlocal] = buf[m++]; - vest[nlocal][0] = buf[m++]; - vest[nlocal][1] = buf[m++]; - vest[nlocal][2] = buf[m++]; - - double **extra = atom->extra; - if (atom->nextra_store) { - int size = static_cast (buf[0]) - m; - for (int i = 0; i < size; i++) - extra[nlocal][i] = buf[m++]; - } - - atom->nlocal++; - return m; -} - -/* ---------------------------------------------------------------------- - create one atom of itype at coord - set other values to defaults - ------------------------------------------------------------------------- */ - -void AtomVecMDPD::create_atom(int itype, double *coord) { - int nlocal = atom->nlocal; - if (nlocal == nmax) - grow(0); - - tag[nlocal] = 0; - type[nlocal] = itype; - x[nlocal][0] = coord[0]; - x[nlocal][1] = coord[1]; - x[nlocal][2] = coord[2]; - mask[nlocal] = 1; - image[nlocal] = ((imageint) IMGMAX << IMG2BITS) | - ((imageint) IMGMAX << IMGBITS) | IMGMAX; - v[nlocal][0] = 0.0; - v[nlocal][1] = 0.0; - v[nlocal][2] = 0.0; - rho[nlocal] = 0.0; - vest[nlocal][0] = 0.0; - vest[nlocal][1] = 0.0; - vest[nlocal][2] = 0.0; - drho[nlocal] = 0.0; - - atom->nlocal++; -} - -/* ---------------------------------------------------------------------- - unpack one line from Atoms section of data file - initialize other atom quantities - ------------------------------------------------------------------------- */ - -void AtomVecMDPD::data_atom(double *coord, imageint imagetmp, char **values) { - int nlocal = atom->nlocal; - if (nlocal == nmax) grow(0); - - tag[nlocal] = utils::tnumeric(FLERR,values[0],true,lmp); - type[nlocal] = utils::inumeric(FLERR,values[1],true,lmp); - if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes) - error->one(FLERR,"Invalid atom type in Atoms section of data file"); - - x[nlocal][0] = coord[0]; - x[nlocal][1] = coord[1]; - x[nlocal][2] = coord[2]; - - image[nlocal] = imagetmp; - - mask[nlocal] = 1; - v[nlocal][0] = 0.0; - v[nlocal][1] = 0.0; - v[nlocal][2] = 0.0; - - vest[nlocal][0] = 0.0; - vest[nlocal][1] = 0.0; - vest[nlocal][2] = 0.0; - - rho[nlocal] = utils::numeric(FLERR,values[2],true,lmp); - drho[nlocal] = 0.0; - - atom->nlocal++; -} - -/* ---------------------------------------------------------------------- - unpack hybrid quantities from one line in Atoms section of data file - initialize other atom quantities for this sub-style - ------------------------------------------------------------------------- */ - -int AtomVecMDPD::data_atom_hybrid(int nlocal, char **values) -{ - rho[nlocal] = utils::numeric(FLERR,values[0],true,lmp); - return 3; -} - -/* ---------------------------------------------------------------------- - pack atom info for data file including 3 image flags + modify what AtomVec::data_atom() just unpacked + or initialize other atom quantities ------------------------------------------------------------------------- */ -void AtomVecMDPD::pack_data(double **buf) +void AtomVecMDPD::data_atom_post(int ilocal) { - int nlocal = atom->nlocal; - for (int i = 0; i < nlocal; i++) { - buf[i][0] = ubuf(tag[i]).d; - buf[i][1] = ubuf(type[i]).d; - buf[i][2] = rho[i]; - buf[i][3] = x[i][0]; - buf[i][4] = x[i][1]; - buf[i][5] = x[i][2]; - buf[i][6] = ubuf((image[i] & IMGMASK) - IMGMAX).d; - buf[i][7] = ubuf((image[i] >> IMGBITS & IMGMASK) - IMGMAX).d; - buf[i][8] = ubuf((image[i] >> IMG2BITS) - IMGMAX).d; - } -} - -/* ---------------------------------------------------------------------- - pack hybrid atom info for data file -------------------------------------------------------------------------- */ - -int AtomVecMDPD::pack_data_hybrid(int i, double *buf) -{ - buf[0] = rho[i]; - return 3; -} - -/* ---------------------------------------------------------------------- - write atom info to data file including 3 image flags -------------------------------------------------------------------------- */ - -void AtomVecMDPD::write_data(FILE *fp, int n, double **buf) -{ - for (int i = 0; i < n; i++) - fprintf(fp,TAGINT_FORMAT - " %d %-1.16e %-1.16e %-1.16e %-1.16e " - "%d %d %d\n", - (tagint) ubuf(buf[i][0]).i,(int) ubuf(buf[i][1]).i, - buf[i][2],buf[i][3],buf[i][4],buf[i][5], - (int) ubuf(buf[i][6]).i,(int) ubuf(buf[i][7]).i, - (int) ubuf(buf[i][8]).i); -} - -/* ---------------------------------------------------------------------- - write hybrid atom info to data file -------------------------------------------------------------------------- */ - -int AtomVecMDPD::write_data_hybrid(FILE *fp, double *buf) -{ - fprintf(fp," %-1.16e",buf[0]); - return 3; + atom->drho[ilocal] = 0.0; + atom->vest[ilocal][0] = 0.0; + atom->vest[ilocal][1] = 0.0; + atom->vest[ilocal][2] = 0.0; } /* ---------------------------------------------------------------------- @@ -912,15 +104,17 @@ void AtomVecMDPD::pack_property_atom(int index, double *buf, { int *mask = atom->mask; int nlocal = atom->nlocal; - int n = 0; + int n = 0; if (index == 0) { + double *rho = atom->rho; for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) buf[n] = rho[i]; else buf[n] = 0.0; n += nvalues; } } else if (index == 1) { + double *drho = atom->drho; for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) buf[n] = drho[i]; else buf[n] = 0.0; @@ -928,24 +122,3 @@ void AtomVecMDPD::pack_property_atom(int index, double *buf, } } } - -/* ---------------------------------------------------------------------- - return # of bytes of allocated memory - ------------------------------------------------------------------------- */ - -bigint AtomVecMDPD::memory_usage() { - bigint bytes = 0; - - if (atom->memcheck("tag")) bytes += memory->usage(tag, nmax); - if (atom->memcheck("type")) bytes += memory->usage(type, nmax); - if (atom->memcheck("mask")) bytes += memory->usage(mask, nmax); - if (atom->memcheck("image")) bytes += memory->usage(image, nmax); - if (atom->memcheck("x")) bytes += memory->usage(x, nmax, 3); - if (atom->memcheck("v")) bytes += memory->usage(v, nmax, 3); - if (atom->memcheck("f")) bytes += memory->usage(f, nmax*comm->nthreads, 3); - if (atom->memcheck("rho")) bytes += memory->usage(rho, nmax); - if (atom->memcheck("drho")) bytes += memory->usage(drho, nmax*comm->nthreads); - if (atom->memcheck("vest")) bytes += memory->usage(vest, nmax, 3); - - return bytes; -} diff --git a/src/USER-MESO/atom_vec_mdpd.h b/src/USER-MESO/atom_vec_mdpd.h index 9e9ffcdcf2..0eb4fff2df 100644 --- a/src/USER-MESO/atom_vec_mdpd.h +++ b/src/USER-MESO/atom_vec_mdpd.h @@ -27,50 +27,11 @@ namespace LAMMPS_NS { class AtomVecMDPD : public AtomVec { public: AtomVecMDPD(class LAMMPS *); - ~AtomVecMDPD() {} - void grow(int); - void grow_reset(); - void copy(int, int, int); + void init(); void force_clear(int, size_t); - int pack_comm(int, int *, double *, int, int *); - int pack_comm_vel(int, int *, double *, int, int *); - void unpack_comm(int, int, double *); - void unpack_comm_vel(int, int, double *); - int pack_reverse(int, int, double *); - void unpack_reverse(int, int *, double *); - int pack_comm_hybrid(int, int *, double *); - int unpack_comm_hybrid(int, int, double *); - int pack_border_hybrid(int, int *, double *); - int unpack_border_hybrid(int, int, double *); - int pack_reverse_hybrid(int, int, double *); - int unpack_reverse_hybrid(int, int *, double *); - int pack_border(int, int *, double *, int, int *); - int pack_border_vel(int, int *, double *, int, int *); - void unpack_border(int, int, double *); - void unpack_border_vel(int, int, double *); - int pack_exchange(int, double *); - int unpack_exchange(double *); - int size_restart(); - int pack_restart(int, double *); - int unpack_restart(double *); - void create_atom(int, double *); - void data_atom(double *, imageint, char **); - int data_atom_hybrid(int, char **); - void pack_data(double **); - int pack_data_hybrid(int, double *); - void write_data(FILE *, int, double **); - int write_data_hybrid(FILE *, double *); + void data_atom_post(int); int property_atom(char *); void pack_property_atom(int, double *, int, int); - bigint memory_usage(); - - private: - tagint *tag; - int *type,*mask; - imageint *image; - double **x,**v,**f; - double *rho, *drho; - double **vest; // estimated velocity during force computation }; } diff --git a/src/USER-MESO/atom_vec_tdpd.cpp b/src/USER-MESO/atom_vec_tdpd.cpp index 74ac47066b..55284f69a2 100644 --- a/src/USER-MESO/atom_vec_tdpd.cpp +++ b/src/USER-MESO/atom_vec_tdpd.cpp @@ -14,12 +14,7 @@ #include "atom_vec_tdpd.h" #include #include "atom.h" -#include "comm.h" -#include "domain.h" -#include "modify.h" -#include "fix.h" #include "update.h" -#include "memory.h" #include "error.h" #include "utils.h" @@ -29,29 +24,30 @@ using namespace LAMMPS_NS; AtomVecTDPD::AtomVecTDPD(LAMMPS *lmp) : AtomVec(lmp) { - if(strcmp(update->unit_style,"lj") != 0) - error->all(FLERR,"Atom style edpd requires lj units"); - molecular = 0; mass_type = 1; forceclearflag = 1; - comm_x_only = comm_f_only = 0; - comm->ghost_velocity = 1; - - cc_species = 0; // for now, reset in process_args() - - size_forward = 3 + cc_species + 3; //vest[3] - size_reverse = 3 + cc_species; - size_border = 6 + cc_species + 3; //vest[3] - size_velocity = 3; - // for data_atom, we read id + type + xyz[3] + cc[i] where i=1,cc_species - size_data_atom = 5 + cc_species; - size_data_vel = 4; - xcol_data = 3; - atom->tdpd_flag = 1; atom->vest_flag = 1; + + // strings with peratom variables to include in each AtomVec method + // strings cannot contain fields in corresponding AtomVec default strings + // order of fields in a string does not matter + // except: fields_data_atom & fields_data_vel must match data file + + fields_grow = (char *) "cc cc_flux vest"; + fields_copy = (char *) "cc vest"; + fields_comm = (char *) "cc vest"; + fields_comm_vel = (char *) "cc vest"; + fields_reverse = (char *) "cc_flux"; + fields_border = (char *) "cc vest"; + fields_border_vel = (char *) "cc vest"; + fields_exchange = (char *) "cc vest"; + fields_restart = (char * ) "cc vest"; + fields_create = (char *) "cc vest"; + fields_data_atom = (char *) "id type x cc"; + fields_data_vel = (char *) "id v"; } /* ---------------------------------------------------------------------- @@ -66,812 +62,42 @@ void AtomVecTDPD::process_args(int narg, char **arg) atom->cc_species = utils::inumeric(FLERR,arg[0],false,lmp); cc_species = atom->cc_species; - // reset sizes that depend on cc_species + atom->add_peratom_change_columns("cc",cc_species); + atom->add_peratom_change_columns("cc_species",cc_species); - size_forward = 3 + cc_species + 3; - size_reverse = 3 + cc_species; - size_border = 6 + cc_species + 3; - size_data_atom = 5 + cc_species; + // delay setting up of fields until now + + setup_fields(); +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecTDPD::init() +{ + AtomVec::init(); + + if (strcmp(update->unit_style,"lj") != 0) + error->all(FLERR,"Atom style tdpd requires lj units"); } /* ---------------------------------------------------------------------- - grow atom arrays - n = 0 grows arrays by a chunk - n > 0 allocates arrays to size n + clear extra forces starting at atom N + nbytes = # of bytes to clear for a per-atom vector ------------------------------------------------------------------------- */ -void AtomVecTDPD::grow(int n) -{ - if (n == 0) grow_nmax(); - else nmax = n; - atom->nmax = nmax; - if (nmax < 0 || nmax > MAXSMALLINT) - error->one(FLERR,"Per-processor system is too big"); - - tag = memory->grow(atom->tag,nmax,"atom:tag"); - type = memory->grow(atom->type,nmax,"atom:type"); - mask = memory->grow(atom->mask,nmax,"atom:mask"); - image = memory->grow(atom->image,nmax,"atom:image"); - x = memory->grow(atom->x,nmax,3,"atom:x"); - v = memory->grow(atom->v,nmax,3,"atom:v"); - f = memory->grow(atom->f,nmax*comm->nthreads,3,"atom:f"); - cc = memory->grow(atom->cc,nmax*comm->nthreads,cc_species,"atom:cc"); - cc_flux = memory->grow(atom->cc_flux,nmax*comm->nthreads,cc_species, - "atom:cc_flux"); - vest = memory->grow(atom->vest, nmax, 3, "atom:vest"); - - if (atom->nextra_grow) - for (int iextra = 0; iextra < atom->nextra_grow; iextra++) - modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax); -} - -/* ---------------------------------------------------------------------- - reset local array ptrs -------------------------------------------------------------------------- */ - -void AtomVecTDPD::grow_reset() -{ - tag = atom->tag; type = atom->type; - mask = atom->mask; image = atom->image; - x = atom->x; v = atom->v; f = atom->f; - cc = atom->cc; cc_flux = atom->cc_flux; - vest = atom->vest; -} - -/* ---------------------------------------------------------------------- - copy atom I info to atom J -------------------------------------------------------------------------- */ - -void AtomVecTDPD::copy(int i, int j, int delflag) -{ - tag[j] = tag[i]; - type[j] = type[i]; - mask[j] = mask[i]; - image[j] = image[i]; - x[j][0] = x[i][0]; - x[j][1] = x[i][1]; - x[j][2] = x[i][2]; - v[j][0] = v[i][0]; - v[j][1] = v[i][1]; - v[j][2] = v[i][2]; - - for(int k = 0; k < cc_species; k++) - cc[j][k] = cc[i][k]; - - vest[j][0] = vest[i][0]; - vest[j][1] = vest[i][1]; - vest[j][2] = vest[i][2]; - - if (atom->nextra_grow) - for (int iextra = 0; iextra < atom->nextra_grow; iextra++) - modify->fix[atom->extra_grow[iextra]]->copy_arrays(i,j,delflag); -} - - void AtomVecTDPD::force_clear(int n, size_t nbytes) { - memset(&cc_flux[n][0],0,cc_species*nbytes); -} - - -/* ---------------------------------------------------------------------- */ - -int AtomVecTDPD::pack_comm(int n, int *list, double *buf, - int pbc_flag, int *pbc) -{ - int i,j,k,m; - double dx,dy,dz; - - m = 0; - if (pbc_flag == 0) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0]; - buf[m++] = x[j][1]; - buf[m++] = x[j][2]; - - for(k = 0; k < cc_species; k++) - buf[m++] = cc[j][k]; - - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } else { - if (domain->triclinic == 0) { - dx = pbc[0]*domain->xprd; - dy = pbc[1]*domain->yprd; - dz = pbc[2]*domain->zprd; - } else { - dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz; - dy = pbc[1]*domain->yprd + pbc[3]*domain->yz; - dz = pbc[2]*domain->zprd; - } - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0] + dx; - buf[m++] = x[j][1] + dy; - buf[m++] = x[j][2] + dz; - - for(k = 0; k < cc_species; k++) - buf[m++] = cc[j][k]; - - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecTDPD::pack_comm_vel(int n, int *list, double *buf, - int pbc_flag, int *pbc) -{ - int i,j,k,m; - double dx,dy,dz,dvx,dvy,dvz; - - m = 0; - if (pbc_flag == 0) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0]; - buf[m++] = x[j][1]; - buf[m++] = x[j][2]; - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; - - for(k = 0; k < cc_species; k++) - buf[m++] = cc[j][k]; - - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } else { - if (domain->triclinic == 0) { - dx = pbc[0]*domain->xprd; - dy = pbc[1]*domain->yprd; - dz = pbc[2]*domain->zprd; - } else { - dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz; - dy = pbc[1]*domain->yprd + pbc[3]*domain->yz; - dz = pbc[2]*domain->zprd; - } - if (!deform_vremap) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0] + dx; - buf[m++] = x[j][1] + dy; - buf[m++] = x[j][2] + dz; - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; - - for(k = 0; k < cc_species; k++) - buf[m++] = cc[j][k]; - - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } else { - dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4]; - dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3]; - dvz = pbc[2]*h_rate[2]; - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0] + dx; - buf[m++] = x[j][1] + dy; - buf[m++] = x[j][2] + dz; - if (mask[i] & deform_groupbit) { - buf[m++] = v[j][0] + dvx; - buf[m++] = v[j][1] + dvy; - buf[m++] = v[j][2] + dvz; - } else { - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; - } - for(k = 0; k < cc_species; k++) - buf[m++] = cc[j][k]; - - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecTDPD::unpack_comm(int n, int first, double *buf) -{ - int i,m,last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - x[i][0] = buf[m++]; - x[i][1] = buf[m++]; - x[i][2] = buf[m++]; - - for(int k = 0; k < cc_species; k++) - cc[i][k] = buf[m++]; - - vest[i][0] = buf[m++]; - vest[i][1] = buf[m++]; - vest[i][2] = buf[m++]; - } -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecTDPD::unpack_comm_vel(int n, int first, double *buf) -{ - int i,m,last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - x[i][0] = buf[m++]; - x[i][1] = buf[m++]; - x[i][2] = buf[m++]; - v[i][0] = buf[m++]; - v[i][1] = buf[m++]; - v[i][2] = buf[m++]; - for(int k = 0; k < cc_species; k++) - cc[i][k] = buf[m++]; - - vest[i][0] = buf[m++]; - vest[i][1] = buf[m++]; - vest[i][2] = buf[m++]; - } -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecTDPD::pack_reverse(int n, int first, double *buf) -{ - int i,m,last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - buf[m++] = f[i][0]; - buf[m++] = f[i][1]; - buf[m++] = f[i][2]; - for(int k = 0; k < cc_species; k++) - buf[m++] = cc_flux[i][k]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecTDPD::unpack_reverse(int n, int *list, double *buf) -{ - int i,j,m; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - f[j][0] += buf[m++]; - f[j][1] += buf[m++]; - f[j][2] += buf[m++]; - for(int k = 0; k < cc_species; k++) - cc_flux[j][k] += buf[m++]; - } -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecTDPD::pack_border(int n, int *list, double *buf, - int pbc_flag, int *pbc) -{ - int i,j,m; - double dx,dy,dz; - - m = 0; - if (pbc_flag == 0) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0]; - buf[m++] = x[j][1]; - buf[m++] = x[j][2]; - buf[m++] = ubuf(tag[j]).d; - buf[m++] = ubuf(type[j]).d; - buf[m++] = ubuf(mask[j]).d; - for(int k = 0; k < cc_species; k++) - buf[m++] = cc[j][k]; - - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } else { - if (domain->triclinic == 0) { - dx = pbc[0]*domain->xprd; - dy = pbc[1]*domain->yprd; - dz = pbc[2]*domain->zprd; - } else { - dx = pbc[0]; - dy = pbc[1]; - dz = pbc[2]; - } - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0] + dx; - buf[m++] = x[j][1] + dy; - buf[m++] = x[j][2] + dz; - buf[m++] = ubuf(tag[j]).d; - buf[m++] = ubuf(type[j]).d; - buf[m++] = ubuf(mask[j]).d; - for(int k = 0; k < cc_species; k++) - buf[m++] = cc[j][k]; - - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } - - if (atom->nextra_border) - for (int iextra = 0; iextra < atom->nextra_border; iextra++) - m += modify->fix[atom->extra_border[iextra]]->pack_border(n,list,&buf[m]); - - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecTDPD::pack_border_vel(int n, int *list, double *buf, - int pbc_flag, int *pbc) -{ - int i,j,m; - double dx,dy,dz,dvx,dvy,dvz; - - m = 0; - if (pbc_flag == 0) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0]; - buf[m++] = x[j][1]; - buf[m++] = x[j][2]; - buf[m++] = ubuf(tag[j]).d; - buf[m++] = ubuf(type[j]).d; - buf[m++] = ubuf(mask[j]).d; - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; - for(int k = 0; k < cc_species; k++) - buf[m++] = cc[j][k]; - - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } else { - if (domain->triclinic == 0) { - dx = pbc[0]*domain->xprd; - dy = pbc[1]*domain->yprd; - dz = pbc[2]*domain->zprd; - } else { - dx = pbc[0]; - dy = pbc[1]; - dz = pbc[2]; - } - if (!deform_vremap) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0] + dx; - buf[m++] = x[j][1] + dy; - buf[m++] = x[j][2] + dz; - buf[m++] = ubuf(tag[j]).d; - buf[m++] = ubuf(type[j]).d; - buf[m++] = ubuf(mask[j]).d; - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; - for(int k = 0; k < cc_species; k++) - buf[m++] = cc[j][k]; - - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } else { - dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4]; - dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3]; - dvz = pbc[2]*h_rate[2]; - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0] + dx; - buf[m++] = x[j][1] + dy; - buf[m++] = x[j][2] + dz; - buf[m++] = ubuf(tag[j]).d; - buf[m++] = ubuf(type[j]).d; - buf[m++] = ubuf(mask[j]).d; - if (mask[i] & deform_groupbit) { - buf[m++] = v[j][0] + dvx; - buf[m++] = v[j][1] + dvy; - buf[m++] = v[j][2] + dvz; - } else { - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; - } - for(int k = 0; k < cc_species; k++) - buf[m++] = cc[j][k]; - - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } - } - - if (atom->nextra_border) - for (int iextra = 0; iextra < atom->nextra_border; iextra++) - m += modify->fix[atom->extra_border[iextra]]->pack_border(n,list,&buf[m]); - - return m; -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecTDPD::unpack_border(int n, int first, double *buf) -{ - int i,m,last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - if (i == nmax) grow(0); - x[i][0] = buf[m++]; - x[i][1] = buf[m++]; - x[i][2] = buf[m++]; - tag[i] = (tagint) ubuf(buf[m++]).i; - type[i] = (int) ubuf(buf[m++]).i; - mask[i] = (int) ubuf(buf[m++]).i; - for(int k = 0; k < cc_species; k++) - cc[i][k] = buf[m++]; - - vest[i][0] = buf[m++]; - vest[i][1] = buf[m++]; - vest[i][2] = buf[m++]; - } - - if (atom->nextra_border) - for (int iextra = 0; iextra < atom->nextra_border; iextra++) - m += modify->fix[atom->extra_border[iextra]]-> - unpack_border(n,first,&buf[m]); -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecTDPD::unpack_border_vel(int n, int first, double *buf) -{ - int i,m,last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - if (i == nmax) grow(0); - x[i][0] = buf[m++]; - x[i][1] = buf[m++]; - x[i][2] = buf[m++]; - tag[i] = (tagint) ubuf(buf[m++]).i; - type[i] = (int) ubuf(buf[m++]).i; - mask[i] = (int) ubuf(buf[m++]).i; - v[i][0] = buf[m++]; - v[i][1] = buf[m++]; - v[i][2] = buf[m++]; - for(int k = 0; k < cc_species; k++) - cc[i][k] = buf[m++]; - - vest[i][0] = buf[m++]; - vest[i][1] = buf[m++]; - vest[i][2] = buf[m++]; - } - - if (atom->nextra_border) - for (int iextra = 0; iextra < atom->nextra_border; iextra++) - m += modify->fix[atom->extra_border[iextra]]-> - unpack_border(n,first,&buf[m]); + memset(&atom->cc_flux[n][0],0,cc_species*nbytes); } /* ---------------------------------------------------------------------- - pack data for atom I for sending to another proc - xyz must be 1st 3 values, so comm::exchange() can test on them + modify what AtomVec::data_atom() just unpacked + or initialize other atom quantities ------------------------------------------------------------------------- */ -int AtomVecTDPD::pack_exchange(int i, double *buf) +void AtomVecTDPD::data_atom_post(int ilocal) { - int m = 1; - buf[m++] = x[i][0]; - buf[m++] = x[i][1]; - buf[m++] = x[i][2]; - buf[m++] = v[i][0]; - buf[m++] = v[i][1]; - buf[m++] = v[i][2]; - buf[m++] = ubuf(tag[i]).d; - buf[m++] = ubuf(type[i]).d; - buf[m++] = ubuf(mask[i]).d; - buf[m++] = ubuf(image[i]).d; - for(int k = 0; k < cc_species; k++) - buf[m++] = cc[i][k]; - - buf[m++] = vest[i][0]; - buf[m++] = vest[i][1]; - buf[m++] = vest[i][2]; - - if (atom->nextra_grow) - for (int iextra = 0; iextra < atom->nextra_grow; iextra++) - m += modify->fix[atom->extra_grow[iextra]]->pack_exchange(i,&buf[m]); - - buf[0] = m; - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecTDPD::unpack_exchange(double *buf) -{ - int nlocal = atom->nlocal; - if (nlocal == nmax) grow(0); - - int m = 1; - x[nlocal][0] = buf[m++]; - x[nlocal][1] = buf[m++]; - x[nlocal][2] = buf[m++]; - v[nlocal][0] = buf[m++]; - v[nlocal][1] = buf[m++]; - v[nlocal][2] = buf[m++]; - tag[nlocal] = (tagint) ubuf(buf[m++]).i; - type[nlocal] = (int) ubuf(buf[m++]).i; - mask[nlocal] = (int) ubuf(buf[m++]).i; - image[nlocal] = (imageint) ubuf(buf[m++]).i; - for(int k = 0; k < cc_species; k++) - cc[nlocal][k] = buf[m++]; - - vest[nlocal][0] = buf[m++]; - vest[nlocal][1] = buf[m++]; - vest[nlocal][2] = buf[m++]; - - if (atom->nextra_grow) - for (int iextra = 0; iextra < atom->nextra_grow; iextra++) - m += modify->fix[atom->extra_grow[iextra]]-> - unpack_exchange(nlocal,&buf[m]); - - atom->nlocal++; - return m; -} - -/* ---------------------------------------------------------------------- - size of restart data for all atoms owned by this proc - include extra data stored by fixes -------------------------------------------------------------------------- */ - -int AtomVecTDPD::size_restart() -{ - int i; - - int nlocal = atom->nlocal; - int n = (11 + cc_species + 3) * nlocal; // 11 + cc[i] + vest[3] - - if (atom->nextra_restart) - for (int iextra = 0; iextra < atom->nextra_restart; iextra++) - for (i = 0; i < nlocal; i++) - n += modify->fix[atom->extra_restart[iextra]]->size_restart(i); - - return n; -} - -/* ---------------------------------------------------------------------- - pack atom I's data for restart file including extra quantities - xyz must be 1st 3 values, so that read_restart can test on them - molecular types may be negative, but write as positive -------------------------------------------------------------------------- */ - -int AtomVecTDPD::pack_restart(int i, double *buf) -{ - int m = 1; - buf[m++] = x[i][0]; - buf[m++] = x[i][1]; - buf[m++] = x[i][2]; - buf[m++] = ubuf(tag[i]).d; - buf[m++] = ubuf(type[i]).d; - buf[m++] = ubuf(mask[i]).d; - buf[m++] = ubuf(image[i]).d; - buf[m++] = v[i][0]; - buf[m++] = v[i][1]; - buf[m++] = v[i][2]; - for(int k = 0; k < cc_species; k++) - buf[m++] = cc[i][k]; - - buf[m++] = vest[i][0]; - buf[m++] = vest[i][1]; - buf[m++] = vest[i][2]; - - if (atom->nextra_restart) - for (int iextra = 0; iextra < atom->nextra_restart; iextra++) - m += modify->fix[atom->extra_restart[iextra]]->pack_restart(i,&buf[m]); - - buf[0] = m; - return m; -} - -/* ---------------------------------------------------------------------- - unpack data for one atom from restart file including extra quantities -------------------------------------------------------------------------- */ - -int AtomVecTDPD::unpack_restart(double *buf) -{ - int nlocal = atom->nlocal; - if (nlocal == nmax) { - grow(0); - if (atom->nextra_store) - memory->grow(atom->extra,nmax,atom->nextra_store,"atom:extra"); - } - - int m = 1; - x[nlocal][0] = buf[m++]; - x[nlocal][1] = buf[m++]; - x[nlocal][2] = buf[m++]; - tag[nlocal] = (tagint) ubuf(buf[m++]).i; - type[nlocal] = (int) ubuf(buf[m++]).i; - mask[nlocal] = (int) ubuf(buf[m++]).i; - image[nlocal] = (imageint) ubuf(buf[m++]).i; - v[nlocal][0] = buf[m++]; - v[nlocal][1] = buf[m++]; - v[nlocal][2] = buf[m++]; - for(int k = 0; k < cc_species; k++) - cc[nlocal][k] = buf[m++]; - - vest[nlocal][0] = buf[m++]; - vest[nlocal][1] = buf[m++]; - vest[nlocal][2] = buf[m++]; - - double **extra = atom->extra; - if (atom->nextra_store) { - int size = static_cast (buf[0]) - m; - for (int i = 0; i < size; i++) extra[nlocal][i] = buf[m++]; - } - - atom->nlocal++; - return m; -} - -/* ---------------------------------------------------------------------- - create one atom of itype at coord - set other values to defaults -------------------------------------------------------------------------- */ - -void AtomVecTDPD::create_atom(int itype, double *coord) -{ - int nlocal = atom->nlocal; - if (nlocal == nmax) grow(0); - - tag[nlocal] = 0; - type[nlocal] = itype; - x[nlocal][0] = coord[0]; - x[nlocal][1] = coord[1]; - x[nlocal][2] = coord[2]; - mask[nlocal] = 1; - image[nlocal] = ((imageint) IMGMAX << IMG2BITS) | - ((imageint) IMGMAX << IMGBITS) | IMGMAX; - v[nlocal][0] = 0.0; - v[nlocal][1] = 0.0; - v[nlocal][2] = 0.0; - - for(int k = 0; k < cc_species; k++) - cc[nlocal][k] = 0.0; - - vest[nlocal][0] = 0.0; - vest[nlocal][1] = 0.0; - vest[nlocal][2] = 0.0; - - atom->nlocal++; -} - -/* ---------------------------------------------------------------------- - unpack one line from Atoms section of data file - initialize other atom quantities -------------------------------------------------------------------------- */ - -void AtomVecTDPD::data_atom(double *coord, imageint imagetmp, char **values) -{ - int nlocal = atom->nlocal; - if (nlocal == nmax) grow(0); - - tag[nlocal] = utils::tnumeric(FLERR,values[0],true,lmp); - type[nlocal] = utils::inumeric(FLERR,values[1],true,lmp); - if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes) - error->one(FLERR,"Invalid atom type in Atoms section of data file"); - - x[nlocal][0] = coord[0]; - x[nlocal][1] = coord[1]; - x[nlocal][2] = coord[2]; - - for(int k = 0; k < cc_species; k++) - cc[nlocal][k] = utils::numeric(FLERR,values[5+k],true,lmp); - - image[nlocal] = imagetmp; - - mask[nlocal] = 1; - v[nlocal][0] = 0.0; - v[nlocal][1] = 0.0; - v[nlocal][2] = 0.0; - - vest[nlocal][0] = 0.0; - vest[nlocal][1] = 0.0; - vest[nlocal][2] = 0.0; - - atom->nlocal++; -} - -/* ---------------------------------------------------------------------- - pack atom info for data file including 3 image flags -------------------------------------------------------------------------- */ - -void AtomVecTDPD::pack_data(double **buf) -{ - int nlocal = atom->nlocal; - for (int i = 0; i < nlocal; i++) { - buf[i][0] = ubuf(tag[i]).d; - buf[i][1] = ubuf(type[i]).d; - buf[i][2] = x[i][0]; - buf[i][3] = x[i][1]; - buf[i][4] = x[i][2]; - buf[i][5] = ubuf((image[i] & IMGMASK) - IMGMAX).d; - buf[i][6] = ubuf((image[i] >> IMGBITS & IMGMASK) - IMGMAX).d; - buf[i][7] = ubuf((image[i] >> IMG2BITS) - IMGMAX).d; - for(int k = 0; k < cc_species; k++) - buf[i][8+k] = cc[i][k]; - } -} - -/* ---------------------------------------------------------------------- - write atom info to data file including 3 image flags -------------------------------------------------------------------------- */ - -void AtomVecTDPD::write_data(FILE *fp, int n, double **buf) -{ - for (int i = 0; i < n; i++){ - fprintf(fp,TAGINT_FORMAT " %d %-1.16e %-1.16e %-1.16e %d %d %d", - (tagint) ubuf(buf[i][0]).i,(int) ubuf(buf[i][1]).i, - buf[i][2],buf[i][3],buf[i][4], - (int) ubuf(buf[i][5]).i,(int) ubuf(buf[i][6]).i, - (int) ubuf(buf[i][7]).i); - for(int k = 0; k < cc_species; k++) - fprintf(fp," %-1.16e",buf[i][8+k]); - fprintf(fp,"\n"); - } -} - -/* ---------------------------------------------------------------------- - return # of bytes of allocated memory -------------------------------------------------------------------------- */ - -bigint AtomVecTDPD::memory_usage() -{ - bigint bytes = 0; - - if (atom->memcheck("tag")) bytes += memory->usage(tag,nmax); - if (atom->memcheck("type")) bytes += memory->usage(type,nmax); - if (atom->memcheck("mask")) bytes += memory->usage(mask,nmax); - if (atom->memcheck("image")) bytes += memory->usage(image,nmax); - if (atom->memcheck("x")) bytes += memory->usage(x,nmax,3); - if (atom->memcheck("v")) bytes += memory->usage(v,nmax,3); - if (atom->memcheck("f")) bytes += memory->usage(f,nmax*comm->nthreads,3); - if (atom->memcheck("cc")) bytes += memory->usage(cc,nmax*comm->nthreads,cc_species); - if (atom->memcheck("cc_flux")) bytes += memory->usage(cc_flux,nmax*comm->nthreads,cc_species); - if (atom->memcheck("vest")) bytes += memory->usage(vest, nmax); - - return bytes; + atom->vest[ilocal][0] = 0.0; + atom->vest[ilocal][1] = 0.0; + atom->vest[ilocal][2] = 0.0; } diff --git a/src/USER-MESO/atom_vec_tdpd.h b/src/USER-MESO/atom_vec_tdpd.h index 86e9ae4bb8..7321859fb4 100644 --- a/src/USER-MESO/atom_vec_tdpd.h +++ b/src/USER-MESO/atom_vec_tdpd.h @@ -27,40 +27,12 @@ namespace LAMMPS_NS { class AtomVecTDPD : public AtomVec { public: AtomVecTDPD(class LAMMPS *); - virtual ~AtomVecTDPD() {} void process_args(int, char **); - void grow(int); - void grow_reset(); - void copy(int, int, int); + void init(); void force_clear(int, size_t); - virtual int pack_comm(int, int *, double *, int, int *); - virtual int pack_comm_vel(int, int *, double *, int, int *); - virtual void unpack_comm(int, int, double *); - virtual void unpack_comm_vel(int, int, double *); - int pack_reverse(int, int, double *); - void unpack_reverse(int, int *, double *); - virtual int pack_border(int, int *, double *, int, int *); - virtual int pack_border_vel(int, int *, double *, int, int *); - virtual void unpack_border(int, int, double *); - virtual void unpack_border_vel(int, int, double *); - virtual int pack_exchange(int, double *); - virtual int unpack_exchange(double *); - int size_restart(); - int pack_restart(int, double *); - int unpack_restart(double *); - void create_atom(int, double *); - void data_atom(double *, imageint, char **); - void pack_data(double **); - void write_data(FILE *, int, double **); - bigint memory_usage(); + void data_atom_post(int); protected: - tagint *tag; - int *type,*mask; - imageint *image; - double **x,**v,**f; - double **vest; // store intermediate velocity for using mvv integrator - double **cc,**cc_flux; int cc_species; }; diff --git a/src/USER-SMD/atom_vec_smd.cpp b/src/USER-SMD/atom_vec_smd.cpp index 604504c5a7..0fada3dc2a 100644 --- a/src/USER-SMD/atom_vec_smd.cpp +++ b/src/USER-SMD/atom_vec_smd.cpp @@ -25,13 +25,7 @@ #include "atom_vec_smd.h" #include #include "atom.h" -#include "comm.h" -#include "domain.h" -#include "modify.h" -#include "fix.h" -#include "memory.h" #include "error.h" -#include "utils.h" using namespace LAMMPS_NS; @@ -40,1247 +34,139 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -AtomVecSMD::AtomVecSMD(LAMMPS *lmp) : - AtomVec(lmp) { - molecular = 0; +AtomVecSMD::AtomVecSMD(LAMMPS *lmp) : AtomVec(lmp) +{ + molecular = 0; + mass_type = 1; + forceclearflag = 1; - comm_x_only = 0; - comm_f_only = 0; - size_forward = 6; // variables that are changed by time integration - size_reverse = 4; // f[3] + de - size_border = 31; - size_velocity = 6; // v + vest - size_data_atom = 13; // 7 + 3 x0 + 3 x - size_data_vel = 4; - xcol_data = 11; + atom->smd_flag = 1; - atom->radius_flag = 1; - atom->rmass_flag = 1; - atom->vfrac_flag = 1; - atom->contact_radius_flag = 1; - atom->molecule_flag = 1; - atom->smd_data_9_flag = 1; - atom->e_flag = 1; - atom->vest_flag = 1; - atom->smd_stress_flag = 1; - atom->eff_plastic_strain_flag = 1; - atom->x0_flag = 1; - atom->damage_flag = 1; - atom->eff_plastic_strain_rate_flag = 1; + atom->radius_flag = 1; + atom->rmass_flag = 1; + atom->vfrac_flag = 1; + atom->contact_radius_flag = 1; + atom->molecule_flag = 1; + atom->smd_data_9_flag = 1; + atom->e_flag = 1; + atom->vest_flag = 1; + atom->smd_stress_flag = 1; + atom->eff_plastic_strain_flag = 1; + atom->x0_flag = 1; + atom->damage_flag = 1; + atom->eff_plastic_strain_rate_flag = 1; - forceclearflag = 1; + // strings with peratom variables to include in each AtomVec method + // strings cannot contain fields in corresponding AtomVec default strings + // order of fields in a string does not matter + // except: fields_data_atom & fields_data_vel must match data file - atom->smd_flag = 1; -} + fields_grow = (char *) + "de vfrac rmass x0 radius contact_radius molecule " + "smd_data_9 e vest tlsph_stress " + "eff_plastic_strain eff_plastic_strain_rate damage"; + fields_copy = (char *) + "vfrac rmass x0 radius contact_radius molecule e " + "eff_plastic_strain eff_plastic_strain_rate vest " + "smd_data_9 smd_stress damage"; + fields_comm = (char *) "radius vfrac vest e"; + fields_comm_vel = (char *) "radius vfrac vest e"; + fields_reverse = (char *) "de"; + fields_border = (char *) + "x0 molecule radius rmass vfrac contact_radius e " + "eff_plastic_strain smd_data_9 smd_stress"; + fields_border_vel = (char *) + "x0 molecule radius rmass vfrac contact_radius e " + "eff_plastic_strain smd_data_9 smd_stress vest"; + fields_exchange = (char *) + "x0 molecule radius rmass vfrac contact_radius e " + "eff_plastic_strain eff_plastic_strain_rate smd_data_9 smd_stress " + "vest damage"; + fields_restart = (char *) + "x0 molecule radius rmass vfrac contact_radius e " + "eff_plastic_strain eff_plastic_strain_rate smd_data_9 smd_stress " + "vest damage"; + fields_create = (char *) + "x0 vest vfrac rmass radius contact_radius molecule e " + "eff_plastic_strain eff_plastic_strain_rate smd_data_9 smd_stress damage"; + fields_data_atom = (char *) + "id type molecule vfrac rmass radius contact_radius x"; + fields_data_vel = (char *) "id v vest"; -/* ---------------------------------------------------------------------- */ + // set these array sizes based on defines -void AtomVecSMD::init() { - AtomVec::init(); + atom->add_peratom_change_columns("smd_data_9",NMAT_FULL); + atom->add_peratom_change_columns("smd_stress",NMAT_SYMM); - // do nothing here + setup_fields(); } /* ---------------------------------------------------------------------- - grow atom arrays - n = 0 grows arrays by a chunk - n > 0 allocates arrays to size n - ------------------------------------------------------------------------- */ + clear extra forces starting at atom N + nbytes = # of bytes to clear for a per-atom vector + NOTE: does f need to be re-cleared? +------------------------------------------------------------------------- */ -void AtomVecSMD::grow(int n) { - if (n == 0) - grow_nmax(); - else - nmax = n; - atom->nmax = nmax; - if (nmax < 0 || nmax > MAXSMALLINT) - error->one(FLERR, "Per-processor system is too big"); - - //printf("in grow, nmax is now %d\n", nmax); - - tag = memory->grow(atom->tag, nmax, "atom:tag"); - type = memory->grow(atom->type, nmax, "atom:type"); - mask = memory->grow(atom->mask, nmax, "atom:mask"); - image = memory->grow(atom->image, nmax, "atom:image"); - x = memory->grow(atom->x, nmax, 3, "atom:x"); - v = memory->grow(atom->v, nmax, 3, "atom:v"); - - f = memory->grow(atom->f, nmax * comm->nthreads, 3, "atom:f"); - de = memory->grow(atom->de, nmax * comm->nthreads, "atom:de"); - - vfrac = memory->grow(atom->vfrac, nmax, "atom:vfrac"); - rmass = memory->grow(atom->rmass, nmax, "atom:rmass"); - x0 = memory->grow(atom->x0, nmax, 3, "atom:x0"); - radius = memory->grow(atom->radius, nmax, "atom:radius"); - contact_radius = memory->grow(atom->contact_radius, nmax, "atom:contact_radius"); - molecule = memory->grow(atom->molecule, nmax, "atom:molecule"); - smd_data_9 = memory->grow(atom->smd_data_9, nmax, NMAT_FULL, "atom:defgrad_old"); - e = memory->grow(atom->e, nmax, "atom:e"); - vest = memory->grow(atom->vest, nmax, 3, "atom:vest"); - tlsph_stress = memory->grow(atom->smd_stress, nmax, NMAT_SYMM, "atom:tlsph_stress"); - eff_plastic_strain = memory->grow(atom->eff_plastic_strain, nmax, "atom:eff_plastic_strain"); - eff_plastic_strain_rate = memory->grow(atom->eff_plastic_strain_rate, nmax, "atom:eff_plastic_strain_rate"); - damage = memory->grow(atom->damage, nmax, "atom:damage"); - - if (atom->nextra_grow) - for (int iextra = 0; iextra < atom->nextra_grow; iextra++) - modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax); +void AtomVecSMD::force_clear(int n, size_t nbytes) +{ + memset(&atom->de[n],0,nbytes); + memset(&atom->f[n][0],0,3*nbytes); } /* ---------------------------------------------------------------------- - reset local array ptrs - ------------------------------------------------------------------------- */ + initialize non-zero atom quantities +------------------------------------------------------------------------- */ -void AtomVecSMD::grow_reset() { - tag = atom->tag; - type = atom->type; - mask = atom->mask; - image = atom->image; - x = atom->x; - v = atom->v; - f = atom->f; - radius = atom->radius; - rmass = atom->rmass; +void AtomVecSMD::create_atom_post(int ilocal) +{ + atom->x0[ilocal][0] = atom->x[ilocal][0]; + atom->x0[ilocal][1] = atom->x[ilocal][1]; + atom->x0[ilocal][2] = atom->x[ilocal][2]; - vfrac = atom->vfrac; - x0 = atom->x0; - contact_radius = atom->contact_radius; - molecule = atom->molecule; - smd_data_9 = atom->smd_data_9; - e = atom->e; - de = atom->de; - tlsph_stress = atom->smd_stress; - eff_plastic_strain = atom->eff_plastic_strain; - eff_plastic_strain_rate = atom->eff_plastic_strain_rate; - damage = atom->damage; - vest = atom->vest; + atom->vfrac[ilocal] = 1.0; + atom->rmass[ilocal] = 1.0; + atom->radius[ilocal] = 0.5; + atom->contact_radius[ilocal] = 0.5; + atom->molecule[ilocal] = 1; + + atom->smd_data_9[ilocal][0] = 1.0; // xx + atom->smd_data_9[ilocal][4] = 1.0; // yy + atom->smd_data_9[ilocal][8] = 1.0; // zz } /* ---------------------------------------------------------------------- - copy atom I info to atom J - ------------------------------------------------------------------------- */ + modify what AtomVec::data_atom() just unpacked + or initialize other atom quantities +------------------------------------------------------------------------- */ -void AtomVecSMD::copy(int i, int j, int delflag) { - tag[j] = tag[i]; - type[j] = type[i]; - mask[j] = mask[i]; - image[j] = image[i]; - x[j][0] = x[i][0]; - x[j][1] = x[i][1]; - x[j][2] = x[i][2]; - v[j][0] = v[i][0]; - v[j][1] = v[i][1]; - v[j][2] = v[i][2]; +void AtomVecSMD::data_atom_post(int ilocal) +{ + atom->e[ilocal] = 0.0; + atom->x0[ilocal][0] = atom->x[ilocal][0]; + atom->x0[ilocal][1] = atom->x[ilocal][1]; + atom->x0[ilocal][2] = atom->x[ilocal][2]; - vfrac[j] = vfrac[i]; - rmass[j] = rmass[i]; - x0[j][0] = x0[i][0]; - x0[j][1] = x0[i][1]; - x0[j][2] = x0[i][2]; - radius[j] = radius[i]; - contact_radius[j] = contact_radius[i]; - molecule[j] = molecule[i]; - e[j] = e[i]; - eff_plastic_strain[j] = eff_plastic_strain[i]; - eff_plastic_strain_rate[j] = eff_plastic_strain_rate[i]; - vest[j][0] = vest[i][0]; - vest[j][1] = vest[i][1]; - vest[j][2] = vest[i][2]; + atom->vest[ilocal][0] = 0.0; + atom->vest[ilocal][1] = 0.0; + atom->vest[ilocal][2] = 0.0; - for (int k = 0; k < NMAT_FULL; k++) { - smd_data_9[j][k] = smd_data_9[i][k]; - } + atom->damage[ilocal] = 0.0; - for (int k = 0; k < NMAT_SYMM; k++) { - tlsph_stress[j][k] = tlsph_stress[i][k]; - } + atom->eff_plastic_strain[ilocal] = 0.0; + atom->eff_plastic_strain_rate[ilocal] = 0.0; - damage[j] = damage[i]; + for (int k = 0; k < NMAT_FULL; k++) + atom->smd_data_9[ilocal][k] = 0.0; - if (atom->nextra_grow) - for (int iextra = 0; iextra < atom->nextra_grow; iextra++) - modify->fix[atom->extra_grow[iextra]]->copy_arrays(i, j, delflag); + for (int k = 0; k < NMAT_SYMM; k++) + atom->smd_stress[ilocal][k] = 0.0; + + atom->smd_data_9[ilocal][0] = 1.0; // xx + atom->smd_data_9[ilocal][4] = 1.0; // yy + atom->smd_data_9[ilocal][8] = 1.0; // zz } -/* ---------------------------------------------------------------------- */ -int AtomVecSMD::pack_comm(int /*n*/, int * /*list*/, double * /*buf*/, int /*pbc_flag*/, int * /*pbc*/) { - error->one(FLERR, "atom vec tlsph can only be used with ghost velocities turned on"); - return -1; -} -/* ---------------------------------------------------------------------- */ -int AtomVecSMD::pack_comm_vel(int n, int *list, double *buf, int pbc_flag, int *pbc) { - // communicate quantities to ghosts, which are changed by time-integration AND are required on ghost atoms. - //no need to pack stress or defgrad information here, as these quantities are not required for ghost atoms. - // Inside pair_style tlsph, these quantities are computed and communicated to ghosts. - - // no need to communicate x0 here, as it is not changed by time integration - // if x0 is changed when the ref config is updated, this communication is performed in the fix_integrate/tlsph - // similarily, rmass could be removed here. - // radius should be communicated here for future time-integration of the radius with ulsph (not implemented yet) - int i, j, m; - double dx, dy, dz, dvx, dvy, dvz; - - m = 0; - if (pbc_flag == 0) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0]; - buf[m++] = x[j][1]; - buf[m++] = x[j][2]; //3 - buf[m++] = radius[j]; - buf[m++] = vfrac[j]; // 5 - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; // 8 - - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; // 11 - buf[m++] = e[j]; // 12 - - } - } else { - if (domain->triclinic == 0) { - dx = pbc[0] * domain->xprd; - dy = pbc[1] * domain->yprd; - dz = pbc[2] * domain->zprd; - } else { - dx = pbc[0] * domain->xprd + pbc[5] * domain->xy + pbc[4] * domain->xz; - dy = pbc[1] * domain->yprd + pbc[3] * domain->yz; - dz = pbc[2] * domain->zprd; - } - if (!deform_vremap) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0] + dx; - buf[m++] = x[j][1] + dy; - buf[m++] = x[j][2] + dz; - buf[m++] = radius[j]; - buf[m++] = vfrac[j]; - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; // 8 - - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; // 11 - buf[m++] = e[j]; // 12 - - } - } else { - dvx = pbc[0] * h_rate[0] + pbc[5] * h_rate[5] + pbc[4] * h_rate[4]; - dvy = pbc[1] * h_rate[1] + pbc[3] * h_rate[3]; - dvz = pbc[2] * h_rate[2]; -// printf("\ndvx = %f, dvy=%f, dvz=%f\n", dvx, dvy, dvz); -// printf("dx = %f, dy=%f, dz=%f\n", dx, dy, dz); - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0] + dx; - buf[m++] = x[j][1] + dy; - buf[m++] = x[j][2] + dz; - buf[m++] = radius[j]; - buf[m++] = vfrac[j]; - if (mask[i] & deform_groupbit) { - buf[m++] = v[j][0] + dvx; - buf[m++] = v[j][1] + dvy; - buf[m++] = v[j][2] + dvz; - buf[m++] = vest[j][0] + dvx; - buf[m++] = vest[j][1] + dvy; - buf[m++] = vest[j][2] + dvz; - } else { - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; // 8 - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; // 11 - } - - buf[m++] = e[j]; // 12 - - } - } - } - - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecSMD::pack_comm_hybrid(int n, int *list, double *buf) { - int i, j, m; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = radius[j]; - buf[m++] = vfrac[j]; - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - buf[m++] = e[j]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecSMD::unpack_comm(int /*n*/, int /*first*/, double * /*buf*/) { - error->one(FLERR, "atom vec tlsph can only be used with ghost velocities turned on"); -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecSMD::unpack_comm_vel(int n, int first, double *buf) { - int i, m, last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - x[i][0] = buf[m++]; - x[i][1] = buf[m++]; - x[i][2] = buf[m++]; //3 - radius[i] = buf[m++]; - vfrac[i] = buf[m++]; // 5 - v[i][0] = buf[m++]; - v[i][1] = buf[m++]; - v[i][2] = buf[m++]; // 8 - - vest[i][0] = buf[m++]; - vest[i][1] = buf[m++]; - vest[i][2] = buf[m++]; // 11 - e[i] = buf[m++]; - - } -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecSMD::unpack_comm_hybrid(int n, int first, double *buf) { - int i, m, last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - radius[i] = buf[m++]; - vfrac[i] = buf[m++]; - vest[i][0] = buf[m++]; - vest[i][1] = buf[m++]; - vest[i][2] = buf[m++]; - e[i] = buf[m++]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecSMD::pack_reverse(int n, int first, double *buf) { - int i, m, last; - - printf("in pack_reverse\n"); - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - buf[m++] = f[i][0]; - buf[m++] = f[i][1]; - buf[m++] = f[i][2]; - buf[m++] = de[i]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecSMD::pack_reverse_hybrid(int n, int first, double *buf) { - int i, m, last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - buf[m++] = de[i]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecSMD::unpack_reverse(int n, int *list, double *buf) { - int i, j, m; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - f[j][0] += buf[m++]; - f[j][1] += buf[m++]; - f[j][2] += buf[m++]; - de[j] += buf[m++]; - } -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecSMD::unpack_reverse_hybrid(int n, int *list, double *buf) { - int i, j, m; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - de[j] += buf[m++]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecSMD::pack_border(int /*n*/, int * /*list*/, double * /*buf*/, int /*pbc_flag*/, int * /*pbc*/) { - error->one(FLERR, "atom vec tlsph can only be used with ghost velocities turned on"); - return -1; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecSMD::pack_border_vel(int n, int *list, double *buf, int pbc_flag, int *pbc) { - int i, j, m; - double dx, dy, dz, dvx, dvy, dvz; - - //printf("AtomVecSMD::pack_border_vel\n"); - - m = 0; - if (pbc_flag == 0) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0]; - buf[m++] = x[j][1]; - buf[m++] = x[j][2]; // 3 - buf[m++] = x0[j][0]; - buf[m++] = x0[j][1]; - buf[m++] = x0[j][2]; // 6 - buf[m++] = ubuf(tag[j]).d; - buf[m++] = ubuf(type[j]).d; - buf[m++] = ubuf(mask[j]).d; - buf[m++] = ubuf(molecule[j]).d; // 10 - buf[m++] = radius[j]; - buf[m++] = rmass[j]; - buf[m++] = vfrac[j]; - buf[m++] = contact_radius[j]; - buf[m++] = e[j]; - buf[m++] = eff_plastic_strain[j]; // 16 - - for (int k = 0; k < NMAT_FULL; k++) { - buf[m++] = smd_data_9[j][k]; - } // 25 - - for (int k = 0; k < NMAT_SYMM; k++) { - buf[m++] = tlsph_stress[j][k]; - } // 31 - - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; // 34 - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; // 37 - } - } else { - - if (domain->triclinic == 0) { - dx = pbc[0] * domain->xprd; - dy = pbc[1] * domain->yprd; - dz = pbc[2] * domain->zprd; - } else { - dx = pbc[0]; - dy = pbc[1]; - dz = pbc[2]; - } - if (!deform_vremap) { - //printf("dx = %f\n", dx); - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0] + dx; - buf[m++] = x[j][1] + dy; - buf[m++] = x[j][2] + dz; // 3 - buf[m++] = x0[j][0]; // this is correct - buf[m++] = x0[j][1]; - buf[m++] = x0[j][2]; // 6 - buf[m++] = ubuf(tag[j]).d; - buf[m++] = ubuf(type[j]).d; - buf[m++] = ubuf(mask[j]).d; - buf[m++] = ubuf(molecule[j]).d; // 10 - buf[m++] = radius[j]; - buf[m++] = rmass[j]; - buf[m++] = vfrac[j]; - buf[m++] = contact_radius[j]; - buf[m++] = e[j]; - buf[m++] = eff_plastic_strain[j]; // 17 - - for (int k = 0; k < NMAT_FULL; k++) { - buf[m++] = smd_data_9[j][k]; - } // 26 - - for (int k = 0; k < NMAT_SYMM; k++) { - buf[m++] = tlsph_stress[j][k]; - } // 32 - - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; // 35 - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; // 38 - - } - } else { - dvx = pbc[0] * h_rate[0] + pbc[5] * h_rate[5] + pbc[4] * h_rate[4]; - dvy = pbc[1] * h_rate[1] + pbc[3] * h_rate[3]; - dvz = pbc[2] * h_rate[2]; -// printf("\ndvx = %f, dvy=%f, dvz=%f\n", dvx, dvy, dvz); -// printf("dx = %f, dy=%f, dz=%f\n", dx, dy, dz); - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0] + dx; - buf[m++] = x[j][1] + dy; - buf[m++] = x[j][2] + dz; // 3 - buf[m++] = x0[j][0]; - buf[m++] = x0[j][1]; - buf[m++] = x0[j][2]; // 6 - buf[m++] = ubuf(tag[j]).d; - buf[m++] = ubuf(type[j]).d; - buf[m++] = ubuf(mask[j]).d; - buf[m++] = ubuf(molecule[j]).d; // 10 - buf[m++] = radius[j]; - buf[m++] = rmass[j]; - buf[m++] = vfrac[j]; - buf[m++] = contact_radius[j]; - buf[m++] = e[j]; - buf[m++] = eff_plastic_strain[j]; // 16 - - for (int k = 0; k < NMAT_FULL; k++) { - buf[m++] = smd_data_9[j][k]; - } // 25 - - for (int k = 0; k < NMAT_SYMM; k++) { - buf[m++] = tlsph_stress[j][k]; - } // 31 - - if (mask[i] & deform_groupbit) { - buf[m++] = v[j][0] + dvx; - buf[m++] = v[j][1] + dvy; - buf[m++] = v[j][2] + dvz; // 34 - buf[m++] = vest[j][0] + dvx; - buf[m++] = vest[j][1] + dvy; - buf[m++] = vest[j][2] + dvz; // 37 - - } else { - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; // 34 - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; // 37 - } - - } - } - } - - if (atom->nextra_border) - for (int iextra = 0; iextra < atom->nextra_border; iextra++) - m += modify->fix[atom->extra_border[iextra]]->pack_border(n, list, &buf[m]); - - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecSMD::pack_border_hybrid(int n, int *list, double *buf) { - int i, j, m; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - - buf[m++] = x0[j][0]; - buf[m++] = x0[j][1]; - buf[m++] = x0[j][2]; // 3 - buf[m++] = ubuf(molecule[j]).d; // 4 - buf[m++] = radius[j]; - buf[m++] = rmass[j]; - buf[m++] = vfrac[j]; - buf[m++] = contact_radius[j]; - buf[m++] = e[j]; - buf[m++] = eff_plastic_strain[j]; // 11 - - for (int k = 0; k < NMAT_FULL; k++) { - buf[m++] = smd_data_9[j][k]; - } // 20 - - for (int k = 0; k < NMAT_SYMM; k++) { - buf[m++] = tlsph_stress[j][k]; - } // 26 - - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecSMD::unpack_border(int /*n*/, int /*first*/, double * /*buf*/) { - error->one(FLERR, "atom vec tlsph can only be used with ghost velocities turned on"); -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecSMD::unpack_border_vel(int n, int first, double *buf) { - int i, m, last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - if (i == nmax) - grow(0); - x[i][0] = buf[m++]; - x[i][1] = buf[m++]; - x[i][2] = buf[m++]; // 3 - x0[i][0] = buf[m++]; - x0[i][1] = buf[m++]; - x0[i][2] = buf[m++]; // 6 - tag[i] = (tagint) ubuf(buf[m++]).i; - type[i] = (int) ubuf(buf[m++]).i; - mask[i] = (int) ubuf(buf[m++]).i; - molecule[i] = (tagint) ubuf(buf[m++]).i; // 10 - - radius[i] = buf[m++]; - rmass[i] = buf[m++]; - vfrac[i] = buf[m++]; - contact_radius[i] = buf[m++]; - e[i] = buf[m++]; - eff_plastic_strain[i] = buf[m++]; // 16 - - for (int k = 0; k < NMAT_FULL; k++) { - smd_data_9[i][k] = buf[m++]; - } // 25 - - for (int k = 0; k < NMAT_SYMM; k++) { - tlsph_stress[i][k] = buf[m++]; - } // 31 - - v[i][0] = buf[m++]; - v[i][1] = buf[m++]; - v[i][2] = buf[m++]; // 34 - vest[i][0] = buf[m++]; - vest[i][1] = buf[m++]; - vest[i][2] = buf[m++]; // 37 - } - - if (atom->nextra_border) - for (int iextra = 0; iextra < atom->nextra_border; iextra++) - m += modify->fix[atom->extra_border[iextra]]->unpack_border(n, first, &buf[m]); -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecSMD::unpack_border_hybrid(int n, int first, double *buf) { - int i, m, last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - x0[i][0] = buf[m++]; - x0[i][1] = buf[m++]; - x0[i][2] = buf[m++]; // 3 - molecule[i] = (tagint) ubuf(buf[m++]).i; // 4 - radius[i] = buf[m++]; - rmass[i] = buf[m++]; - vfrac[i] = buf[m++]; - contact_radius[i] = buf[m++]; - e[i] = buf[m++]; - eff_plastic_strain[i] = buf[m++]; // 11 - - for (int k = 0; k < NMAT_FULL; k++) { - smd_data_9[i][k] = buf[m++]; - } // 20 - - for (int k = 0; k < NMAT_SYMM; k++) { - tlsph_stress[i][k] = buf[m++]; - } // 26 - } - return m; -} - -/* ---------------------------------------------------------------------- - pack data for atom I for sending to another proc - xyz must be 1st 3 values, so comm::exchange() can test on them - ------------------------------------------------------------------------- */ - -int AtomVecSMD::pack_exchange(int i, double *buf) { - int m = 1; - - //printf("in AtomVecSMD::pack_exchange tag %d\n", tag[i]); - - buf[m++] = x[i][0]; - buf[m++] = x[i][1]; - buf[m++] = x[i][2]; // 3 - buf[m++] = x0[i][0]; - buf[m++] = x0[i][1]; - buf[m++] = x0[i][2]; // 6 - buf[m++] = ubuf(tag[i]).d; - buf[m++] = ubuf(type[i]).d; - buf[m++] = ubuf(mask[i]).d; - buf[m++] = ubuf(image[i]).d; - buf[m++] = ubuf(molecule[i]).d; // 11 - buf[m++] = radius[i]; - buf[m++] = rmass[i]; - buf[m++] = vfrac[i]; - buf[m++] = contact_radius[i]; - buf[m++] = e[i]; - buf[m++] = eff_plastic_strain[i]; // 18 - buf[m++] = eff_plastic_strain_rate[i]; // 19 - - for (int k = 0; k < NMAT_FULL; k++) { - buf[m++] = smd_data_9[i][k]; - } // 27 - - for (int k = 0; k < NMAT_SYMM; k++) { - buf[m++] = tlsph_stress[i][k]; - } // 33 - - buf[m++] = v[i][0]; - buf[m++] = v[i][1]; - buf[m++] = v[i][2]; // 36 - buf[m++] = vest[i][0]; - buf[m++] = vest[i][1]; - buf[m++] = vest[i][2]; // 39 - - buf[m++] = damage[i]; - - if (atom->nextra_grow) - for (int iextra = 0; iextra < atom->nextra_grow; iextra++) - m += modify->fix[atom->extra_grow[iextra]]->pack_exchange(i, &buf[m]); - - buf[0] = m; - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecSMD::unpack_exchange(double *buf) { - int nlocal = atom->nlocal; - if (nlocal == nmax) - grow(0); - - int m = 1; - - x[nlocal][0] = buf[m++]; - x[nlocal][1] = buf[m++]; - x[nlocal][2] = buf[m++]; // 3 - x0[nlocal][0] = buf[m++]; - x0[nlocal][1] = buf[m++]; - x0[nlocal][2] = buf[m++]; // 6 - tag[nlocal] = (tagint) ubuf(buf[m++]).i; - type[nlocal] = (int) ubuf(buf[m++]).i; - mask[nlocal] = (int) ubuf(buf[m++]).i; - image[nlocal] = (imageint) ubuf(buf[m++]).i; - molecule[nlocal] = (tagint) ubuf(buf[m++]).i; // 11 - - radius[nlocal] = buf[m++]; - rmass[nlocal] = buf[m++]; - vfrac[nlocal] = buf[m++]; - contact_radius[nlocal] = buf[m++]; - e[nlocal] = buf[m++]; - eff_plastic_strain[nlocal] = buf[m++]; // 18 - eff_plastic_strain_rate[nlocal] = buf[m++]; // 19 - - for (int k = 0; k < NMAT_FULL; k++) { - smd_data_9[nlocal][k] = buf[m++]; - } // 27 - - for (int k = 0; k < NMAT_SYMM; k++) { - tlsph_stress[nlocal][k] = buf[m++]; - } // 33 - - v[nlocal][0] = buf[m++]; - v[nlocal][1] = buf[m++]; - v[nlocal][2] = buf[m++]; // 36 - vest[nlocal][0] = buf[m++]; - vest[nlocal][1] = buf[m++]; - vest[nlocal][2] = buf[m++]; // 39 - - damage[nlocal] = buf[m++]; //40 - - if (atom->nextra_grow) - for (int iextra = 0; iextra < atom->nextra_grow; iextra++) - m += modify->fix[atom->extra_grow[iextra]]->unpack_exchange(nlocal, &buf[m]); - - atom->nlocal++; - return m; -} - -/* ---------------------------------------------------------------------- - size of restart data for all atoms owned by this proc - include extra data stored by fixes - ------------------------------------------------------------------------- */ - -int AtomVecSMD::size_restart() { - int i; - - int nlocal = atom->nlocal; - int n = 43 * nlocal; // count pack_restart + 1 (size of buffer) - - if (atom->nextra_restart) - for (int iextra = 0; iextra < atom->nextra_restart; iextra++) - for (i = 0; i < nlocal; i++) - n += modify->fix[atom->extra_restart[iextra]]->size_restart(i); - - return n; -} - -/* ---------------------------------------------------------------------- - pack atom I's data for restart file including extra quantities - xyz must be 1st 3 values, so that read_restart can test on them - molecular types may be negative, but write as positive - ------------------------------------------------------------------------- */ -int AtomVecSMD::pack_restart(int i, double *buf) { - int m = 1; // 1 - - buf[m++] = x[i][0]; - buf[m++] = x[i][1]; - buf[m++] = x[i][2]; // 4 - buf[m++] = x0[i][0]; - buf[m++] = x0[i][1]; - buf[m++] = x0[i][2]; // 7 - buf[m++] = ubuf(tag[i]).d; - buf[m++] = ubuf(type[i]).d; - buf[m++] = ubuf(mask[i]).d; // 10 - buf[m++] = ubuf(image[i]).d; - buf[m++] = ubuf(molecule[i]).d; - buf[m++] = radius[i]; - buf[m++] = rmass[i]; - buf[m++] = vfrac[i]; // 15 - buf[m++] = contact_radius[i]; - buf[m++] = e[i]; - buf[m++] = eff_plastic_strain[i]; - buf[m++] = eff_plastic_strain_rate[i]; // 19 - - for (int k = 0; k < NMAT_FULL; k++) { - buf[m++] = smd_data_9[i][k]; - } // 28 - - for (int k = 0; k < NMAT_SYMM; k++) { - buf[m++] = tlsph_stress[i][k]; - } // 34 - - buf[m++] = v[i][0]; - buf[m++] = v[i][1]; - buf[m++] = v[i][2]; // 37 - buf[m++] = vest[i][0]; - buf[m++] = vest[i][1]; - buf[m++] = vest[i][2]; // 40 - - buf[m++] = damage[i]; // 41 - - if (atom->nextra_restart) - for (int iextra = 0; iextra < atom->nextra_restart; iextra++) - m += modify->fix[atom->extra_restart[iextra]]->pack_restart(i, &buf[m]); - - buf[0] = m; - return m; -} - -/* ---------------------------------------------------------------------- - unpack data for one atom from restart file including extra quantities - ------------------------------------------------------------------------- */ - -int AtomVecSMD::unpack_restart(double *buf) { - int nlocal = atom->nlocal; - if (nlocal == nmax) { - grow(0); - if (atom->nextra_store) - memory->grow(atom->extra, nmax, atom->nextra_store, "atom:extra"); - } - - int m = 1; - - x[nlocal][0] = buf[m++]; - x[nlocal][1] = buf[m++]; - x[nlocal][2] = buf[m++]; // 3 - x0[nlocal][0] = buf[m++]; - x0[nlocal][1] = buf[m++]; - x0[nlocal][2] = buf[m++]; // 6 - tag[nlocal] = (tagint) ubuf(buf[m++]).i; - type[nlocal] = (int) ubuf(buf[m++]).i; - mask[nlocal] = (int) ubuf(buf[m++]).i; - image[nlocal] = (imageint) ubuf(buf[m++]).i; - molecule[nlocal] = (tagint) ubuf(buf[m++]).i; // 11 - - radius[nlocal] = buf[m++]; - rmass[nlocal] = buf[m++]; - vfrac[nlocal] = buf[m++]; //14 - contact_radius[nlocal] = buf[m++]; //15 - e[nlocal] = buf[m++]; - eff_plastic_strain[nlocal] = buf[m++]; // 18 - eff_plastic_strain_rate[nlocal] = buf[m++]; // 29 - - for (int k = 0; k < NMAT_FULL; k++) { - smd_data_9[nlocal][k] = buf[m++]; - } // 28 - - for (int k = 0; k < NMAT_SYMM; k++) { - tlsph_stress[nlocal][k] = buf[m++]; - } // 34 - - v[nlocal][0] = buf[m++]; - v[nlocal][1] = buf[m++]; - v[nlocal][2] = buf[m++]; // 37 - vest[nlocal][0] = buf[m++]; - vest[nlocal][1] = buf[m++]; - vest[nlocal][2] = buf[m++]; // 40 - - damage[nlocal] = buf[m++]; //41 - - //printf("nlocal in restart is %d\n", nlocal); - - double **extra = atom->extra; - if (atom->nextra_store) { - int size = static_cast(buf[0]) - m; - for (int i = 0; i < size; i++) - extra[nlocal][i] = buf[m++]; - } - - atom->nlocal++; - - //printf("returning m=%d in unpack_restart\n", m); - - return m; -} - -/* ---------------------------------------------------------------------- - create one atom of itype at coord - set other values to defaults - ------------------------------------------------------------------------- */ - -void AtomVecSMD::create_atom(int itype, double *coord) { - int nlocal = atom->nlocal; - if (nlocal == nmax) { - printf("nlocal = %d, nmax = %d, calling grow\n", nlocal, nmax); - grow(0); - printf("... finished growing\n"); - } - - tag[nlocal] = 0; - type[nlocal] = itype; - x[nlocal][0] = coord[0]; - x[nlocal][1] = coord[1]; - x[nlocal][2] = coord[2]; - x0[nlocal][0] = coord[0]; - x0[nlocal][1] = coord[1]; - x0[nlocal][2] = coord[2]; - mask[nlocal] = 1; - image[nlocal] = ((imageint) IMGMAX << IMG2BITS) | ((imageint) IMGMAX << IMGBITS) | IMGMAX; - v[nlocal][0] = 0.0; - v[nlocal][1] = 0.0; - v[nlocal][2] = 0.0; - vest[nlocal][0] = 0.0; - vest[nlocal][1] = 0.0; - vest[nlocal][2] = 0.0; - - vfrac[nlocal] = 1.0; - rmass[nlocal] = 1.0; - radius[nlocal] = 0.5; - contact_radius[nlocal] = 0.5; - molecule[nlocal] = 1; - e[nlocal] = 0.0; - eff_plastic_strain[nlocal] = 0.0; - eff_plastic_strain_rate[nlocal] = 0.0; - - for (int k = 0; k < NMAT_FULL; k++) { - smd_data_9[nlocal][k] = 0.0; - } - smd_data_9[nlocal][0] = 1.0; // xx - smd_data_9[nlocal][4] = 1.0; // yy - smd_data_9[nlocal][8] = 1.0; // zz - - for (int k = 0; k < NMAT_SYMM; k++) { - tlsph_stress[nlocal][k] = 0.0; - } - - damage[nlocal] = 0.0; - - atom->nlocal++; -} - -/* ---------------------------------------------------------------------- - unpack one line from Atoms section of data file - initialize other atom quantities - ------------------------------------------------------------------------- */ - -void AtomVecSMD::data_atom(double *coord, imageint imagetmp, char **values) { - int nlocal = atom->nlocal; - if (nlocal == nmax) - grow(0); - - tag[nlocal] = utils::tnumeric(FLERR,values[0],true,lmp); - - type[nlocal] = utils::inumeric(FLERR,values[1],true,lmp); - if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes) - error->one(FLERR, "Invalid atom type in Atoms section of data file"); - - molecule[nlocal] = utils::tnumeric(FLERR,values[2],true,lmp); - if (molecule[nlocal] <= 0) - error->one(FLERR, "Invalid molecule in Atoms section of data file"); - - vfrac[nlocal] = utils::numeric(FLERR,values[3],true,lmp); - if (vfrac[nlocal] < 0.0) - error->one(FLERR, "Invalid volume in Atoms section of data file"); - - rmass[nlocal] = utils::numeric(FLERR,values[4],true,lmp); - if (rmass[nlocal] == 0.0) - error->one(FLERR, "Invalid mass in Atoms section of data file"); - - radius[nlocal] = utils::numeric(FLERR,values[5],true,lmp); - if (radius[nlocal] < 0.0) - error->one(FLERR, "Invalid radius in Atoms section of data file"); - - contact_radius[nlocal] = utils::numeric(FLERR,values[6],true,lmp); - if (contact_radius[nlocal] < 0.0) - error->one(FLERR, "Invalid contact radius in Atoms section of data file"); - - e[nlocal] = 0.0; - - x0[nlocal][0] = utils::numeric(FLERR,values[7],true,lmp); - x0[nlocal][1] = utils::numeric(FLERR,values[8],true,lmp); - x0[nlocal][2] = utils::numeric(FLERR,values[9],true,lmp); - - x[nlocal][0] = coord[0]; - x[nlocal][1] = coord[1]; - x[nlocal][2] = coord[2]; - - image[nlocal] = imagetmp; - - mask[nlocal] = 1; - v[nlocal][0] = 0.0; - v[nlocal][1] = 0.0; - v[nlocal][2] = 0.0; - vest[nlocal][0] = 0.0; - vest[nlocal][1] = 0.0; - vest[nlocal][2] = 0.0; - - damage[nlocal] = 0.0; - - eff_plastic_strain[nlocal] = 0.0; - eff_plastic_strain_rate[nlocal] = 0.0; - - for (int k = 0; k < NMAT_FULL; k++) { - smd_data_9[nlocal][k] = 0.0; - } - - for (int k = 0; k < NMAT_SYMM; k++) { - tlsph_stress[nlocal][k] = 0.0; - } - - smd_data_9[nlocal][0] = 1.0; // xx - smd_data_9[nlocal][4] = 1.0; // yy - smd_data_9[nlocal][8] = 1.0; // zz - - atom->nlocal++; -} - -/* ---------------------------------------------------------------------- - unpack hybrid quantities from one line in Atoms section of data file - initialize other atom quantities for this sub-style - ------------------------------------------------------------------------- */ - -int AtomVecSMD::data_atom_hybrid(int /*nlocal*/, char **/*values*/) { - error->one(FLERR, "hybrid atom style functionality not yet implemented for atom style smd"); - return -1; -} - -/* ---------------------------------------------------------------------- - unpack one line from Velocities section of data file - ------------------------------------------------------------------------- */ - -void AtomVecSMD::data_vel(int m, char **values) { - v[m][0] = utils::numeric(FLERR,values[0],true,lmp); - v[m][1] = utils::numeric(FLERR,values[1],true,lmp); - v[m][2] = utils::numeric(FLERR,values[2],true,lmp); - vest[m][0] = utils::numeric(FLERR,values[0],true,lmp); - vest[m][1] = utils::numeric(FLERR,values[1],true,lmp); - vest[m][2] = utils::numeric(FLERR,values[2],true,lmp); -} - -/* ---------------------------------------------------------------------- - unpack hybrid quantities from one line in Velocities section of data file - ------------------------------------------------------------------------- */ - -int AtomVecSMD::data_vel_hybrid(int /*m*/, char **/*values*/) { - error->one(FLERR, "hybrid atom style functionality not yet implemented for atom style smd"); - return 0; -} - -/* ---------------------------------------------------------------------- - pack atom info for data file including 3 image flags - ------------------------------------------------------------------------- */ - -void AtomVecSMD::pack_data(double **buf) { - int nlocal = atom->nlocal; - for (int i = 0; i < nlocal; i++) { - buf[i][0] = ubuf(tag[i]).d; - buf[i][1] = ubuf(type[i]).d; - buf[i][2] = ubuf(molecule[i]).d; - buf[i][3] = vfrac[i]; - buf[i][4] = rmass[i]; - buf[i][5] = radius[i]; - buf[i][6] = contact_radius[i]; - - buf[i][7] = x[i][0]; - buf[i][8] = x[i][1]; - buf[i][9] = x[i][2]; - - buf[i][10] = x0[i][0]; - buf[i][11] = x0[i][1]; - buf[i][12] = x0[i][2]; - - buf[i][13] = ubuf((image[i] & IMGMASK) - IMGMAX).d; - buf[i][14] = ubuf((image[i] >> IMGBITS & IMGMASK) - IMGMAX).d; - buf[i][15] = ubuf((image[i] >> IMG2BITS) - IMGMAX).d; - } -} - -/* ---------------------------------------------------------------------- - pack hybrid atom info for data file - ------------------------------------------------------------------------- */ - -int AtomVecSMD::pack_data_hybrid(int /*i*/, double * /*buf*/) { - error->one(FLERR, "hybrid atom style functionality not yet implemented for atom style smd"); - return -1; -} - -/* ---------------------------------------------------------------------- - write atom info to data file including 3 image flags - ------------------------------------------------------------------------- */ - -void AtomVecSMD::write_data(FILE *fp, int n, double **buf) { - for (int i = 0; i < n; i++) - fprintf(fp, - TAGINT_FORMAT - " %d %d %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e\n", (tagint) ubuf(buf[i][0]).i, - (int) ubuf(buf[i][1]).i, (int) ubuf(buf[i][2]).i, buf[i][3], buf[i][4], buf[i][5], buf[i][6], buf[i][7], buf[i][8], - buf[i][9], buf[i][10], buf[i][11], buf[i][12]); -} - -/* ---------------------------------------------------------------------- - write hybrid atom info to data file - ------------------------------------------------------------------------- */ - -int AtomVecSMD::write_data_hybrid(FILE * /*fp*/, double * /*buf*/) { - error->one(FLERR, "hybrid atom style functionality not yet implemented for atom style smd"); - return -1; -} - -/* ---------------------------------------------------------------------- - pack velocity info for data file - ------------------------------------------------------------------------- */ - -void AtomVecSMD::pack_vel(double **buf) { - int nlocal = atom->nlocal; - for (int i = 0; i < nlocal; i++) { - buf[i][0] = ubuf(tag[i]).d; - buf[i][1] = v[i][0]; - buf[i][2] = v[i][1]; - buf[i][3] = v[i][2]; - } -} - -/* ---------------------------------------------------------------------- - pack hybrid velocity info for data file - ------------------------------------------------------------------------- */ - -int AtomVecSMD::pack_vel_hybrid(int /*i*/, double * /*buf*/) { - error->one(FLERR, "hybrid atom style functionality not yet implemented for atom style smd"); - return 0; -} - -/* ---------------------------------------------------------------------- - write velocity info to data file - ------------------------------------------------------------------------- */ - -void AtomVecSMD::write_vel(FILE *fp, int n, double **buf) { - for (int i = 0; i < n; i++) - fprintf(fp, TAGINT_FORMAT - " %-1.16e %-1.16e %-1.16e\n", (tagint) ubuf(buf[i][0]).i, buf[i][1], buf[i][2], buf[i][3]); -} - -/* ---------------------------------------------------------------------- - write hybrid velocity info to data file - ------------------------------------------------------------------------- */ - -int AtomVecSMD::write_vel_hybrid(FILE * /*fp*/, double * /*buf*/) { - error->one(FLERR, "hybrid atom style functionality not yet implemented for atom style smd"); - return 3; -} - -/* ---------------------------------------------------------------------- - return # of bytes of allocated memory - ------------------------------------------------------------------------- */ - -bigint AtomVecSMD::memory_usage() { - bigint bytes = 0; - - if (atom->memcheck("tag")) - bytes += memory->usage(tag, nmax); - if (atom->memcheck("type")) - bytes += memory->usage(type, nmax); - if (atom->memcheck("molecule")) - bytes += memory->usage(molecule, nmax); - if (atom->memcheck("mask")) - bytes += memory->usage(mask, nmax); - if (atom->memcheck("image")) - bytes += memory->usage(image, nmax); - if (atom->memcheck("x")) - bytes += memory->usage(x, nmax, 3); - if (atom->memcheck("v")) - bytes += memory->usage(v, nmax, 3); - if (atom->memcheck("vest")) - bytes += memory->usage(vest, nmax, 3); - if (atom->memcheck("f")) - bytes += memory->usage(f, nmax * comm->nthreads, 3); - - if (atom->memcheck("radius")) - bytes += memory->usage(radius, nmax); - if (atom->memcheck("contact_radius")) - bytes += memory->usage(contact_radius, nmax); - if (atom->memcheck("vfrac")) - bytes += memory->usage(vfrac, nmax); - if (atom->memcheck("rmass")) - bytes += memory->usage(rmass, nmax); - if (atom->memcheck("eff_plastic_strain")) - bytes += memory->usage(eff_plastic_strain, nmax); - if (atom->memcheck("eff_plastic_strain_rate")) - bytes += memory->usage(eff_plastic_strain_rate, nmax); - if (atom->memcheck("e")) - bytes += memory->usage(e, nmax); - if (atom->memcheck("de")) - bytes += memory->usage(de, nmax); - - if (atom->memcheck("smd_data_9")) - bytes += memory->usage(smd_data_9, nmax, NMAT_FULL); - if (atom->memcheck("tlsph_stress")) - bytes += memory->usage(tlsph_stress, nmax, NMAT_SYMM); - - if (atom->memcheck("damage")) - bytes += memory->usage(damage, nmax); - - return bytes; -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecSMD::force_clear(int n, size_t nbytes) { - //printf("clearing force on atom %d", n); - memset(&de[n], 0, nbytes); - memset(&f[0][0], 0, 3 * nbytes); -} diff --git a/src/USER-SMD/atom_vec_smd.h b/src/USER-SMD/atom_vec_smd.h index 34fdfc1f76..539f209ca7 100644 --- a/src/USER-SMD/atom_vec_smd.h +++ b/src/USER-SMD/atom_vec_smd.h @@ -9,7 +9,6 @@ * * ----------------------------------------------------------------------- */ - /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories @@ -39,63 +38,9 @@ namespace LAMMPS_NS { class AtomVecSMD : public AtomVec { public: AtomVecSMD(class LAMMPS *); - ~AtomVecSMD() {} - void init(); - void grow(int); - void grow_reset(); - void copy(int, int, int); void force_clear(int, size_t); - int pack_comm(int, int *, double *, int, int *); - int pack_comm_vel(int, int *, double *, int, int *); - int pack_comm_hybrid(int, int *, double *); - void unpack_comm(int, int, double *); - void unpack_comm_vel(int, int, double *); - int unpack_comm_hybrid(int, int, double *); - int pack_reverse(int, int, double *); - int pack_reverse_hybrid(int, int, double *); - void unpack_reverse(int, int *, double *); - int unpack_reverse_hybrid(int, int *, double *); - int pack_border(int, int *, double *, int, int *); - int pack_border_vel(int, int *, double *, int, int *); - int pack_border_hybrid(int, int *, double *); - void unpack_border(int, int, double *); - void unpack_border_vel(int, int, double *); - int unpack_border_hybrid(int, int, double *); - int pack_exchange(int, double *); - int unpack_exchange(double *); - int size_restart(); - int pack_restart(int, double *); - int unpack_restart(double *); - void create_atom(int, double *); - void data_atom(double *, imageint, char **); - int data_atom_hybrid(int, char **); - void data_vel(int, char **); - int data_vel_hybrid(int, char **); - void pack_data(double **); - int pack_data_hybrid(int, double *); - void write_data(FILE *, int, double **); - int write_data_hybrid(FILE *, double *); - void pack_vel(double **); - int pack_vel_hybrid(int, double *); - void write_vel(FILE *, int, double **); - int write_vel_hybrid(FILE *, double *); - bigint memory_usage(); - - private: - tagint *tag; - int *type,*mask; - imageint *image; - double **x,**v,**f; - double *radius,*rmass; - - tagint *molecule; - double *vfrac,**x0,*contact_radius, **smd_data_9, *e, *de, **vest; - double **tlsph_stress; - double *eff_plastic_strain; - double *damage; - double *eff_plastic_strain_rate; - - + void create_atom_post(int); + void data_atom_post(int); }; } diff --git a/src/USER-SPH/atom_vec_meso.cpp b/src/USER-SPH/atom_vec_meso.cpp index cd7c2251ab..a80ab91d2e 100644 --- a/src/USER-SPH/atom_vec_meso.cpp +++ b/src/USER-SPH/atom_vec_meso.cpp @@ -14,13 +14,7 @@ #include "atom_vec_meso.h" #include #include "atom.h" -#include "comm.h" -#include "domain.h" -#include "modify.h" -#include "fix.h" -#include "memory.h" #include "error.h" -#include "utils.h" using namespace LAMMPS_NS; @@ -32,921 +26,64 @@ AtomVecMeso::AtomVecMeso(LAMMPS *lmp) : AtomVec(lmp) mass_type = 1; forceclearflag = 1; - comm_x_only = 0; // we communicate not only x forward but also vest ... - comm_f_only = 0; // we also communicate de and drho in reverse direction - size_forward = 8; // 3 + rho + e + vest[3], that means we may only communicate 5 in hybrid - size_reverse = 5; // 3 + drho + de - size_border = 12; // 6 + rho + e + vest[3] + cv - size_velocity = 3; - size_data_atom = 8; - size_data_vel = 4; - xcol_data = 6; - atom->e_flag = 1; atom->rho_flag = 1; atom->cv_flag = 1; atom->vest_flag = 1; + + // strings with peratom variables to include in each AtomVec method + // strings cannot contain fields in corresponding AtomVec default strings + // order of fields in a string does not matter + // except: fields_data_atom & fields_data_vel must match data file + + fields_grow = (char *) "rho drho e de cv vest"; + fields_copy = (char *) "rho drho e de cv vest"; + fields_comm = (char *) "rho e vest"; + fields_comm_vel = (char *) "rho e vest"; + fields_reverse = (char *) "drho de"; + fields_border = (char *) "rho e cv vest"; + fields_border_vel = (char *) "rho e cv vest"; + fields_exchange = (char *) "rho e cv vest"; + fields_restart = (char * ) "rho e cv vest"; + fields_create = (char *) "rho e cv vest de drho"; + fields_data_atom = (char *) "id type rho e cv x"; + fields_data_vel = (char *) "id v"; + + setup_fields(); } /* ---------------------------------------------------------------------- - grow atom arrays - n = 0 grows arrays by a chunk - n > 0 allocates arrays to size n - ------------------------------------------------------------------------- */ - -void AtomVecMeso::grow(int n) -{ - if (n == 0) grow_nmax(); - else nmax = n; - atom->nmax = nmax; - if (nmax < 0 || nmax > MAXSMALLINT) - error->one(FLERR,"Per-processor system is too big"); - - tag = memory->grow(atom->tag, nmax, "atom:tag"); - type = memory->grow(atom->type, nmax, "atom:type"); - mask = memory->grow(atom->mask, nmax, "atom:mask"); - image = memory->grow(atom->image, nmax, "atom:image"); - x = memory->grow(atom->x, nmax, 3, "atom:x"); - v = memory->grow(atom->v, nmax, 3, "atom:v"); - f = memory->grow(atom->f, nmax*comm->nthreads, 3, "atom:f"); - - rho = memory->grow(atom->rho, nmax, "atom:rho"); - drho = memory->grow(atom->drho, nmax*comm->nthreads, "atom:drho"); - e = memory->grow(atom->e, nmax, "atom:e"); - de = memory->grow(atom->de, nmax*comm->nthreads, "atom:de"); - vest = memory->grow(atom->vest, nmax, 3, "atom:vest"); - cv = memory->grow(atom->cv, nmax, "atom:cv"); - - if (atom->nextra_grow) - for (int iextra = 0; iextra < atom->nextra_grow; iextra++) - modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax); -} - -/* ---------------------------------------------------------------------- - reset local array ptrs - ------------------------------------------------------------------------- */ - -void AtomVecMeso::grow_reset() { - tag = atom->tag; - type = atom->type; - mask = atom->mask; - image = atom->image; - x = atom->x; - v = atom->v; - f = atom->f; - rho = atom->rho; - drho = atom->drho; - e = atom->e; - de = atom->de; - vest = atom->vest; - cv = atom->cv; -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecMeso::copy(int i, int j, int delflag) { - //printf("in AtomVecMeso::copy\n"); - tag[j] = tag[i]; - type[j] = type[i]; - mask[j] = mask[i]; - image[j] = image[i]; - x[j][0] = x[i][0]; - x[j][1] = x[i][1]; - x[j][2] = x[i][2]; - v[j][0] = v[i][0]; - v[j][1] = v[i][1]; - v[j][2] = v[i][2]; - - rho[j] = rho[i]; - drho[j] = drho[i]; - e[j] = e[i]; - de[j] = de[i]; - cv[j] = cv[i]; - vest[j][0] = vest[i][0]; - vest[j][1] = vest[i][1]; - vest[j][2] = vest[i][2]; - - if (atom->nextra_grow) - for (int iextra = 0; iextra < atom->nextra_grow; iextra++) - modify->fix[atom->extra_grow[iextra]]->copy_arrays(i, j,delflag); -} - -/* ---------------------------------------------------------------------- */ + clear extra forces starting at atom N + nbytes = # of bytes to clear for a per-atom vector +------------------------------------------------------------------------- */ void AtomVecMeso::force_clear(int n, size_t nbytes) { - memset(&de[n],0,nbytes); - memset(&drho[n],0,nbytes); -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecMeso::pack_comm_hybrid(int n, int *list, double *buf) { - //printf("in AtomVecMeso::pack_comm_hybrid\n"); - int i, j, m; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = rho[j]; - buf[m++] = e[j]; - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecMeso::unpack_comm_hybrid(int n, int first, double *buf) { - //printf("in AtomVecMeso::unpack_comm_hybrid\n"); - int i, m, last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - rho[i] = buf[m++]; - e[i] = buf[m++]; - vest[i][0] = buf[m++]; - vest[i][1] = buf[m++]; - vest[i][2] = buf[m++]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecMeso::pack_border_hybrid(int n, int *list, double *buf) { - //printf("in AtomVecMeso::pack_border_hybrid\n"); - int i, j, m; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = rho[j]; - buf[m++] = e[j]; - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecMeso::unpack_border_hybrid(int n, int first, double *buf) { - //printf("in AtomVecMeso::unpack_border_hybrid\n"); - int i, m, last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - rho[i] = buf[m++]; - e[i] = buf[m++]; - vest[i][0] = buf[m++]; - vest[i][1] = buf[m++]; - vest[i][2] = buf[m++]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecMeso::pack_reverse_hybrid(int n, int first, double *buf) { - //printf("in AtomVecMeso::pack_reverse_hybrid\n"); - int i, m, last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - buf[m++] = drho[i]; - buf[m++] = de[i]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecMeso::unpack_reverse_hybrid(int n, int *list, double *buf) { - //printf("in AtomVecMeso::unpack_reverse_hybrid\n"); - int i, j, m; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - drho[j] += buf[m++]; - de[j] += buf[m++]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecMeso::pack_comm(int n, int *list, double *buf, int pbc_flag, - int *pbc) { - //printf("in AtomVecMeso::pack_comm\n"); - int i, j, m; - double dx, dy, dz; - - m = 0; - if (pbc_flag == 0) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0]; - buf[m++] = x[j][1]; - buf[m++] = x[j][2]; - buf[m++] = rho[j]; - buf[m++] = e[j]; - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } else { - if (domain->triclinic == 0) { - dx = pbc[0] * domain->xprd; - dy = pbc[1] * domain->yprd; - dz = pbc[2] * domain->zprd; - } else { - dx = pbc[0] * domain->xprd + pbc[5] * domain->xy + pbc[4] * domain->xz; - dy = pbc[1] * domain->yprd + pbc[3] * domain->yz; - dz = pbc[2] * domain->zprd; - } - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0] + dx; - buf[m++] = x[j][1] + dy; - buf[m++] = x[j][2] + dz; - buf[m++] = rho[j]; - buf[m++] = e[j]; - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecMeso::pack_comm_vel(int n, int *list, double *buf, int pbc_flag, - int *pbc) { - //printf("in AtomVecMeso::pack_comm_vel\n"); - int i, j, m; - double dx, dy, dz; - - m = 0; - if (pbc_flag == 0) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0]; - buf[m++] = x[j][1]; - buf[m++] = x[j][2]; - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; - buf[m++] = rho[j]; - buf[m++] = e[j]; - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } else { - if (domain->triclinic == 0) { - dx = pbc[0] * domain->xprd; - dy = pbc[1] * domain->yprd; - dz = pbc[2] * domain->zprd; - } else { - dx = pbc[0] * domain->xprd + pbc[5] * domain->xy + pbc[4] * domain->xz; - dy = pbc[1] * domain->yprd + pbc[3] * domain->yz; - dz = pbc[2] * domain->zprd; - } - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0] + dx; - buf[m++] = x[j][1] + dy; - buf[m++] = x[j][2] + dz; - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; - buf[m++] = rho[j]; - buf[m++] = e[j]; - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecMeso::unpack_comm(int n, int first, double *buf) { - //printf("in AtomVecMeso::unpack_comm\n"); - int i, m, last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - x[i][0] = buf[m++]; - x[i][1] = buf[m++]; - x[i][2] = buf[m++]; - rho[i] = buf[m++]; - e[i] = buf[m++]; - vest[i][0] = buf[m++]; - vest[i][1] = buf[m++]; - vest[i][2] = buf[m++]; - } -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecMeso::unpack_comm_vel(int n, int first, double *buf) { - //printf("in AtomVecMeso::unpack_comm_vel\n"); - int i, m, last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - x[i][0] = buf[m++]; - x[i][1] = buf[m++]; - x[i][2] = buf[m++]; - v[i][0] = buf[m++]; - v[i][1] = buf[m++]; - v[i][2] = buf[m++]; - rho[i] = buf[m++]; - e[i] = buf[m++]; - vest[i][0] = buf[m++]; - vest[i][1] = buf[m++]; - vest[i][2] = buf[m++]; - } -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecMeso::pack_reverse(int n, int first, double *buf) { - //printf("in AtomVecMeso::pack_reverse\n"); - int i, m, last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - buf[m++] = f[i][0]; - buf[m++] = f[i][1]; - buf[m++] = f[i][2]; - buf[m++] = drho[i]; - buf[m++] = de[i]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecMeso::unpack_reverse(int n, int *list, double *buf) { - //printf("in AtomVecMeso::unpack_reverse\n"); - int i, j, m; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - f[j][0] += buf[m++]; - f[j][1] += buf[m++]; - f[j][2] += buf[m++]; - drho[j] += buf[m++]; - de[j] += buf[m++]; - } -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecMeso::pack_border(int n, int *list, double *buf, int pbc_flag, - int *pbc) { - //printf("in AtomVecMeso::pack_border\n"); - int i, j, m; - double dx, dy, dz; - - m = 0; - if (pbc_flag == 0) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0]; - buf[m++] = x[j][1]; - buf[m++] = x[j][2]; - buf[m++] = ubuf(tag[j]).d; - buf[m++] = ubuf(type[j]).d; - buf[m++] = ubuf(mask[j]).d; - buf[m++] = rho[j]; - buf[m++] = e[j]; - buf[m++] = cv[j]; - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } else { - if (domain->triclinic == 0) { - dx = pbc[0] * domain->xprd; - dy = pbc[1] * domain->yprd; - dz = pbc[2] * domain->zprd; - } else { - dx = pbc[0]; - dy = pbc[1]; - dz = pbc[2]; - } - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0] + dx; - buf[m++] = x[j][1] + dy; - buf[m++] = x[j][2] + dz; - buf[m++] = ubuf(tag[j]).d; - buf[m++] = ubuf(type[j]).d; - buf[m++] = ubuf(mask[j]).d; - buf[m++] = rho[j]; - buf[m++] = e[j]; - buf[m++] = cv[j]; - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } - - if (atom->nextra_border) - for (int iextra = 0; iextra < atom->nextra_border; iextra++) - m += modify->fix[atom->extra_border[iextra]]->pack_border(n,list,&buf[m]); - - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecMeso::pack_border_vel(int n, int *list, double *buf, int pbc_flag, - int *pbc) -{ - int i,j,m; - double dx,dy,dz,dvx,dvy,dvz; - - m = 0; - if (pbc_flag == 0) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0]; - buf[m++] = x[j][1]; - buf[m++] = x[j][2]; - buf[m++] = ubuf(tag[j]).d; - buf[m++] = ubuf(type[j]).d; - buf[m++] = ubuf(mask[j]).d; - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; - buf[m++] = rho[j]; - buf[m++] = e[j]; - buf[m++] = cv[j]; - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } else { - if (domain->triclinic == 0) { - dx = pbc[0] * domain->xprd; - dy = pbc[1] * domain->yprd; - dz = pbc[2] * domain->zprd; - } else { - dx = pbc[0]; - dy = pbc[1]; - dz = pbc[2]; - } - if (!deform_vremap) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0] + dx; - buf[m++] = x[j][1] + dy; - buf[m++] = x[j][2] + dz; - buf[m++] = ubuf(tag[j]).d; - buf[m++] = ubuf(type[j]).d; - buf[m++] = ubuf(mask[j]).d; - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; - buf[m++] = rho[j]; - buf[m++] = e[j]; - buf[m++] = cv[j]; - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - } else { - dvx = pbc[0] * h_rate[0] + pbc[5] * h_rate[5] + pbc[4] * h_rate[4]; - dvy = pbc[1] * h_rate[1] + pbc[3] * h_rate[3]; - dvz = pbc[2] * h_rate[2]; - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = x[j][0] + dx; - buf[m++] = x[j][1] + dy; - buf[m++] = x[j][2] + dz; - buf[m++] = ubuf(tag[j]).d; - buf[m++] = ubuf(type[j]).d; - buf[m++] = ubuf(mask[j]).d; - if (mask[i] & deform_groupbit) { - buf[m++] = v[j][0] + dvx; - buf[m++] = v[j][1] + dvy; - buf[m++] = v[j][2] + dvz; - buf[m++] = vest[j][0] + dvx; - buf[m++] = vest[j][1] + dvy; - buf[m++] = vest[j][2] + dvz; - } else { - buf[m++] = v[j][0]; - buf[m++] = v[j][1]; - buf[m++] = v[j][2]; - buf[m++] = vest[j][0]; - buf[m++] = vest[j][1]; - buf[m++] = vest[j][2]; - } - buf[m++] = rho[j]; - buf[m++] = e[j]; - buf[m++] = cv[j]; - } - } - } - - if (atom->nextra_border) - for (int iextra = 0; iextra < atom->nextra_border; iextra++) - m += modify->fix[atom->extra_border[iextra]]->pack_border(n,list,&buf[m]); - - return m; -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecMeso::unpack_border(int n, int first, double *buf) { - //printf("in AtomVecMeso::unpack_border\n"); - int i, m, last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - if (i == nmax) - grow(0); - x[i][0] = buf[m++]; - x[i][1] = buf[m++]; - x[i][2] = buf[m++]; - tag[i] = (tagint) ubuf(buf[m++]).i; - type[i] = (int) ubuf(buf[m++]).i; - mask[i] = (int) ubuf(buf[m++]).i; - rho[i] = buf[m++]; - e[i] = buf[m++]; - cv[i] = buf[m++]; - vest[i][0] = buf[m++]; - vest[i][1] = buf[m++]; - vest[i][2] = buf[m++]; - } - - if (atom->nextra_border) - for (int iextra = 0; iextra < atom->nextra_border; iextra++) - m += modify->fix[atom->extra_border[iextra]]-> - unpack_border(n,first,&buf[m]); -} - -/* ---------------------------------------------------------------------- */ - -void AtomVecMeso::unpack_border_vel(int n, int first, double *buf) { - //printf("in AtomVecMeso::unpack_border_vel\n"); - int i, m, last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - if (i == nmax) - grow(0); - x[i][0] = buf[m++]; - x[i][1] = buf[m++]; - x[i][2] = buf[m++]; - tag[i] = (tagint) ubuf(buf[m++]).i; - type[i] = (int) ubuf(buf[m++]).i; - mask[i] = (int) ubuf(buf[m++]).i; - v[i][0] = buf[m++]; - v[i][1] = buf[m++]; - v[i][2] = buf[m++]; - vest[i][0] = buf[m++]; - vest[i][1] = buf[m++]; - vest[i][2] = buf[m++]; - rho[i] = buf[m++]; - e[i] = buf[m++]; - cv[i] = buf[m++]; - } - - if (atom->nextra_border) - for (int iextra = 0; iextra < atom->nextra_border; iextra++) - m += modify->fix[atom->extra_border[iextra]]-> - unpack_border(n,first,&buf[m]); + memset(&atom->de[n],0,nbytes); + memset(&atom->drho[n],0,nbytes); } /* ---------------------------------------------------------------------- - pack data for atom I for sending to another proc - xyz must be 1st 3 values, so comm::exchange() can test on them - ------------------------------------------------------------------------- */ - -int AtomVecMeso::pack_exchange(int i, double *buf) { - //printf("in AtomVecMeso::pack_exchange\n"); - int m = 1; - buf[m++] = x[i][0]; - buf[m++] = x[i][1]; - buf[m++] = x[i][2]; - buf[m++] = v[i][0]; - buf[m++] = v[i][1]; - buf[m++] = v[i][2]; - buf[m++] = ubuf(tag[i]).d; - buf[m++] = ubuf(type[i]).d; - buf[m++] = ubuf(mask[i]).d; - buf[m++] = ubuf(image[i]).d; - buf[m++] = rho[i]; - buf[m++] = e[i]; - buf[m++] = cv[i]; - buf[m++] = vest[i][0]; - buf[m++] = vest[i][1]; - buf[m++] = vest[i][2]; - - if (atom->nextra_grow) - for (int iextra = 0; iextra < atom->nextra_grow; iextra++) - m += modify->fix[atom->extra_grow[iextra]]->pack_exchange(i, &buf[m]); - - buf[0] = m; - return m; -} - -/* ---------------------------------------------------------------------- */ - -int AtomVecMeso::unpack_exchange(double *buf) { - //printf("in AtomVecMeso::unpack_exchange\n"); - int nlocal = atom->nlocal; - if (nlocal == nmax) - grow(0); - - int m = 1; - x[nlocal][0] = buf[m++]; - x[nlocal][1] = buf[m++]; - x[nlocal][2] = buf[m++]; - v[nlocal][0] = buf[m++]; - v[nlocal][1] = buf[m++]; - v[nlocal][2] = buf[m++]; - tag[nlocal] = (tagint) ubuf(buf[m++]).i; - type[nlocal] = (int) ubuf(buf[m++]).i; - mask[nlocal] = (int) ubuf(buf[m++]).i; - image[nlocal] = (imageint) ubuf(buf[m++]).i; - rho[nlocal] = buf[m++]; - e[nlocal] = buf[m++]; - cv[nlocal] = buf[m++]; - vest[nlocal][0] = buf[m++]; - vest[nlocal][1] = buf[m++]; - vest[nlocal][2] = buf[m++]; - - if (atom->nextra_grow) - for (int iextra = 0; iextra < atom->nextra_grow; iextra++) - m += modify->fix[atom->extra_grow[iextra]]-> unpack_exchange(nlocal, - &buf[m]); - - atom->nlocal++; - return m; -} - -/* ---------------------------------------------------------------------- - size of restart data for all atoms owned by this proc - include extra data stored by fixes - ------------------------------------------------------------------------- */ - -int AtomVecMeso::size_restart() { - int i; - - int nlocal = atom->nlocal; - int n = 17 * nlocal; // 11 + rho + e + cv + vest[3] - - if (atom->nextra_restart) - for (int iextra = 0; iextra < atom->nextra_restart; iextra++) - for (i = 0; i < nlocal; i++) - n += modify->fix[atom->extra_restart[iextra]]->size_restart(i); - - return n; -} - -/* ---------------------------------------------------------------------- - pack atom I's data for restart file including extra quantities - xyz must be 1st 3 values, so that read_restart can test on them - molecular types may be negative, but write as positive - ------------------------------------------------------------------------- */ - -int AtomVecMeso::pack_restart(int i, double *buf) { - int m = 1; - buf[m++] = x[i][0]; - buf[m++] = x[i][1]; - buf[m++] = x[i][2]; - buf[m++] = ubuf(tag[i]).d; - buf[m++] = ubuf(type[i]).d; - buf[m++] = ubuf(mask[i]).d; - buf[m++] = ubuf(image[i]).d; - buf[m++] = v[i][0]; - buf[m++] = v[i][1]; - buf[m++] = v[i][2]; - buf[m++] = rho[i]; - buf[m++] = e[i]; - buf[m++] = cv[i]; - buf[m++] = vest[i][0]; - buf[m++] = vest[i][1]; - buf[m++] = vest[i][2]; - - if (atom->nextra_restart) - for (int iextra = 0; iextra < atom->nextra_restart; iextra++) - m += modify->fix[atom->extra_restart[iextra]]->pack_restart(i, &buf[m]); - - buf[0] = m; - return m; -} - -/* ---------------------------------------------------------------------- - unpack data for one atom from restart file including extra quantities - ------------------------------------------------------------------------- */ - -int AtomVecMeso::unpack_restart(double *buf) { - int nlocal = atom->nlocal; - if (nlocal == nmax) { - grow(0); - if (atom->nextra_store) - memory->grow(atom->extra, nmax, atom->nextra_store, "atom:extra"); - } - - int m = 1; - x[nlocal][0] = buf[m++]; - x[nlocal][1] = buf[m++]; - x[nlocal][2] = buf[m++]; - tag[nlocal] = (tagint) ubuf(buf[m++]).i; - type[nlocal] = (int) ubuf(buf[m++]).i; - mask[nlocal] = (int) ubuf(buf[m++]).i; - image[nlocal] = (imageint) ubuf(buf[m++]).i; - v[nlocal][0] = buf[m++]; - v[nlocal][1] = buf[m++]; - v[nlocal][2] = buf[m++]; - rho[nlocal] = buf[m++]; - e[nlocal] = buf[m++]; - cv[nlocal] = buf[m++]; - vest[nlocal][0] = buf[m++]; - vest[nlocal][1] = buf[m++]; - vest[nlocal][2] = buf[m++]; - - double **extra = atom->extra; - if (atom->nextra_store) { - int size = static_cast (buf[0]) - m; - for (int i = 0; i < size; i++) - extra[nlocal][i] = buf[m++]; - } - - atom->nlocal++; - return m; -} - -/* ---------------------------------------------------------------------- - create one atom of itype at coord - set other values to defaults - ------------------------------------------------------------------------- */ - -void AtomVecMeso::create_atom(int itype, double *coord) { - int nlocal = atom->nlocal; - if (nlocal == nmax) - grow(0); - - tag[nlocal] = 0; - type[nlocal] = itype; - x[nlocal][0] = coord[0]; - x[nlocal][1] = coord[1]; - x[nlocal][2] = coord[2]; - mask[nlocal] = 1; - image[nlocal] = ((imageint) IMGMAX << IMG2BITS) | - ((imageint) IMGMAX << IMGBITS) | IMGMAX; - v[nlocal][0] = 0.0; - v[nlocal][1] = 0.0; - v[nlocal][2] = 0.0; - rho[nlocal] = 0.0; - e[nlocal] = 0.0; - cv[nlocal] = 1.0; - vest[nlocal][0] = 0.0; - vest[nlocal][1] = 0.0; - vest[nlocal][2] = 0.0; - de[nlocal] = 0.0; - drho[nlocal] = 0.0; - - atom->nlocal++; -} - -/* ---------------------------------------------------------------------- - unpack one line from Atoms section of data file - initialize other atom quantities - ------------------------------------------------------------------------- */ - -void AtomVecMeso::data_atom(double *coord, imageint imagetmp, char **values) { - int nlocal = atom->nlocal; - if (nlocal == nmax) grow(0); - - tag[nlocal] = utils::tnumeric(FLERR,values[0],true,lmp); - type[nlocal] = utils::inumeric(FLERR,values[1],true,lmp); - if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes) - error->one(FLERR,"Invalid atom type in Atoms section of data file"); - - rho[nlocal] = utils::numeric(FLERR,values[2],true,lmp); - e[nlocal] = utils::numeric(FLERR,values[3],true,lmp); - cv[nlocal] = utils::numeric(FLERR,values[4],true,lmp); - - x[nlocal][0] = coord[0]; - x[nlocal][1] = coord[1]; - x[nlocal][2] = coord[2]; - - //printf("rho=%f, e=%f, cv=%f, x=%f\n", rho[nlocal], e[nlocal], cv[nlocal], x[nlocal][0]); - - image[nlocal] = imagetmp; - - mask[nlocal] = 1; - v[nlocal][0] = 0.0; - v[nlocal][1] = 0.0; - v[nlocal][2] = 0.0; - - vest[nlocal][0] = 0.0; - vest[nlocal][1] = 0.0; - vest[nlocal][2] = 0.0; - - de[nlocal] = 0.0; - drho[nlocal] = 0.0; - - atom->nlocal++; -} - -/* ---------------------------------------------------------------------- - unpack hybrid quantities from one line in Atoms section of data file - initialize other atom quantities for this sub-style - ------------------------------------------------------------------------- */ - -int AtomVecMeso::data_atom_hybrid(int nlocal, char **values) { - - rho[nlocal] = utils::numeric(FLERR,values[0],true,lmp); - e[nlocal] = utils::numeric(FLERR,values[1],true,lmp); - cv[nlocal] = utils::numeric(FLERR,values[2],true,lmp); - - return 3; -} - -/* ---------------------------------------------------------------------- - pack atom info for data file including 3 image flags + initialize non-zero atom quantities ------------------------------------------------------------------------- */ -void AtomVecMeso::pack_data(double **buf) +void AtomVecMeso::create_atom_post(int ilocal) { - int nlocal = atom->nlocal; - for (int i = 0; i < nlocal; i++) { - buf[i][0] = ubuf(tag[i]).d; - buf[i][1] = ubuf(type[i]).d; - buf[i][2] = rho[i]; - buf[i][3] = e[i]; - buf[i][4] = cv[i]; - buf[i][5] = x[i][0]; - buf[i][6] = x[i][1]; - buf[i][7] = x[i][2]; - buf[i][8] = ubuf((image[i] & IMGMASK) - IMGMAX).d; - buf[i][9] = ubuf((image[i] >> IMGBITS & IMGMASK) - IMGMAX).d; - buf[i][10] = ubuf((image[i] >> IMG2BITS) - IMGMAX).d; - } + atom->cv[ilocal] = 1.0; } /* ---------------------------------------------------------------------- - pack hybrid atom info for data file + modify what AtomVec::data_atom() just unpacked + or initialize other atom quantities ------------------------------------------------------------------------- */ -int AtomVecMeso::pack_data_hybrid(int i, double *buf) +void AtomVecMeso::data_atom_post(int ilocal) { - buf[0] = rho[i]; - buf[1] = e[i]; - buf[2] = cv[i]; - return 3; -} - -/* ---------------------------------------------------------------------- - write atom info to data file including 3 image flags -------------------------------------------------------------------------- */ - -void AtomVecMeso::write_data(FILE *fp, int n, double **buf) -{ - for (int i = 0; i < n; i++) - fprintf(fp,TAGINT_FORMAT - " %d %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e " - "%d %d %d\n", - (tagint) ubuf(buf[i][0]).i,(int) ubuf(buf[i][1]).i, - buf[i][2],buf[i][3],buf[i][4], - buf[i][5],buf[i][6],buf[i][7], - (int) ubuf(buf[i][8]).i,(int) ubuf(buf[i][9]).i, - (int) ubuf(buf[i][10]).i); -} - -/* ---------------------------------------------------------------------- - write hybrid atom info to data file -------------------------------------------------------------------------- */ - -int AtomVecMeso::write_data_hybrid(FILE *fp, double *buf) -{ - fprintf(fp," %-1.16e %-1.16e %-1.16e",buf[0],buf[1],buf[2]); - return 3; + atom->vest[ilocal][0] = 0.0; + atom->vest[ilocal][1] = 0.0; + atom->vest[ilocal][2] = 0.0; + atom->de[ilocal] = 0.0; + atom->drho[ilocal] = 0.0; } /* ---------------------------------------------------------------------- @@ -977,30 +114,35 @@ void AtomVecMeso::pack_property_atom(int index, double *buf, int n = 0; if (index == 0) { + double *rho = atom->rho; for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) buf[n] = rho[i]; else buf[n] = 0.0; n += nvalues; } } else if (index == 1) { + double *drho = atom->drho; for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) buf[n] = drho[i]; else buf[n] = 0.0; n += nvalues; } } else if (index == 2) { + double *e = atom->e; for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) buf[n] = e[i]; else buf[n] = 0.0; n += nvalues; } } else if (index == 3) { + double *de = atom->de; for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) buf[n] = de[i]; else buf[n] = 0.0; n += nvalues; } } else if (index == 4) { + double *cv = atom->cv; for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) buf[n] = cv[i]; else buf[n] = 0.0; @@ -1008,40 +150,3 @@ void AtomVecMeso::pack_property_atom(int index, double *buf, } } } - -/* ---------------------------------------------------------------------- - return # of bytes of allocated memory - ------------------------------------------------------------------------- */ - -bigint AtomVecMeso::memory_usage() { - bigint bytes = 0; - - if (atom->memcheck("tag")) - bytes += memory->usage(tag, nmax); - if (atom->memcheck("type")) - bytes += memory->usage(type, nmax); - if (atom->memcheck("mask")) - bytes += memory->usage(mask, nmax); - if (atom->memcheck("image")) - bytes += memory->usage(image, nmax); - if (atom->memcheck("x")) - bytes += memory->usage(x, nmax, 3); - if (atom->memcheck("v")) - bytes += memory->usage(v, nmax, 3); - if (atom->memcheck("f")) - bytes += memory->usage(f, nmax*comm->nthreads, 3); - if (atom->memcheck("rho")) - bytes += memory->usage(rho, nmax); - if (atom->memcheck("drho")) - bytes += memory->usage(drho, nmax*comm->nthreads); - if (atom->memcheck("e")) - bytes += memory->usage(e, nmax); - if (atom->memcheck("de")) - bytes += memory->usage(de, nmax*comm->nthreads); - if (atom->memcheck("cv")) - bytes += memory->usage(cv, nmax); - if (atom->memcheck("vest")) - bytes += memory->usage(vest, nmax); - - return bytes; -} diff --git a/src/USER-SPH/atom_vec_meso.h b/src/USER-SPH/atom_vec_meso.h index da68222e29..c8a1090474 100644 --- a/src/USER-SPH/atom_vec_meso.h +++ b/src/USER-SPH/atom_vec_meso.h @@ -27,50 +27,11 @@ namespace LAMMPS_NS { class AtomVecMeso : public AtomVec { public: AtomVecMeso(class LAMMPS *); - ~AtomVecMeso() {} - void grow(int); - void grow_reset(); - void copy(int, int, int); void force_clear(int, size_t); - int pack_comm(int, int *, double *, int, int *); - int pack_comm_vel(int, int *, double *, int, int *); - void unpack_comm(int, int, double *); - void unpack_comm_vel(int, int, double *); - int pack_reverse(int, int, double *); - void unpack_reverse(int, int *, double *); - int pack_comm_hybrid(int, int *, double *); - int unpack_comm_hybrid(int, int, double *); - int pack_border_hybrid(int, int *, double *); - int unpack_border_hybrid(int, int, double *); - int pack_reverse_hybrid(int, int, double *); - int unpack_reverse_hybrid(int, int *, double *); - int pack_border(int, int *, double *, int, int *); - int pack_border_vel(int, int *, double *, int, int *); - void unpack_border(int, int, double *); - void unpack_border_vel(int, int, double *); - int pack_exchange(int, double *); - int unpack_exchange(double *); - int size_restart(); - int pack_restart(int, double *); - int unpack_restart(double *); - void create_atom(int, double *); - void data_atom(double *, imageint, char **); - int data_atom_hybrid(int, char **); - void pack_data(double **); - int pack_data_hybrid(int, double *); - void write_data(FILE *, int, double **); - int write_data_hybrid(FILE *, double *); + void create_atom_post(int); + void data_atom_post(int); int property_atom(char *); void pack_property_atom(int, double *, int, int); - bigint memory_usage(); - - private: - tagint *tag; - int *type,*mask; - imageint *image; - double **x,**v,**f; - double *rho, *drho, *e, *de, *cv; - double **vest; // estimated velocity during force computation }; } diff --git a/src/atom.cpp b/src/atom.cpp index d238db5e65..8fd859478f 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -81,50 +81,21 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) image = NULL; x = v = f = NULL; - molecule = NULL; - molindex = molatom = NULL; + // charged and dipolar particles + q = NULL; mu = NULL; + + // finite-size particles + omega = angmom = torque = NULL; radius = rmass = NULL; ellipsoid = line = tri = body = NULL; - vfrac = s0 = NULL; - x0 = NULL; + // molecular systems - spin = NULL; - eradius = ervel = erforce = NULL; - cs = csforce = vforce = ervelforce = NULL; - etag = NULL; - - rho = drho = e = de = cv = NULL; - vest = NULL; - - // SPIN package - - sp = fm = fm_long = NULL; - - // USER-DPD - - uCond = uMech = uChem = uCG = uCGnew = NULL; - duChem = NULL; - dpdTheta = NULL; - - // USER-MESO - - cc = cc_flux = NULL; - edpd_temp = edpd_flux = edpd_cv = NULL; - - // USER-SMD - - contact_radius = NULL; - smd_data_9 = NULL; - smd_stress = NULL; - eff_plastic_strain = NULL; - eff_plastic_strain_rate = NULL; - damage = NULL; - - // molecular info + molecule = NULL; + molindex = molatom = NULL; bond_per_atom = extra_bond_per_atom = 0; num_bond = NULL; @@ -150,6 +121,47 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) nspecial = NULL; special = NULL; + // PERI package + + vfrac = s0 = NULL; + x0 = NULL; + + // SPIN package + + sp = fm = fm_long = NULL; + + // USER-EFF and USER-AWPMD packages + + spin = NULL; + eradius = ervel = erforce = NULL; + ervelforce = cs = csforce = NULL; + vforce = NULL; + etag = NULL; + + // USER-DPD package + + uCond = uMech = uChem = uCG = uCGnew = NULL; + duChem = dpdTheta = NULL; + + // USER-MESO package + + cc = cc_flux = NULL; + edpd_temp = edpd_flux = edpd_cv = NULL; + + // USER-SMD package + + contact_radius = NULL; + smd_data_9 = NULL; + smd_stress = NULL; + eff_plastic_strain = NULL; + eff_plastic_strain_rate = NULL; + damage = NULL; + + // USER-SPH package + + rho = drho = e = de = cv = NULL; + vest = NULL; + // user-defined molecules nmolecule = 0; @@ -532,9 +544,16 @@ void Atom::peratom_create() add_peratom("ervel",&ervel,DOUBLE,0); add_peratom("erforce",&erforce,DOUBLE,0,1); // set per-thread flag + // USER-AWPMD package + + add_peratom("cs",&cs,DOUBLE,0); + add_peratom("csforce",&csforce,DOUBLE,0); + add_peratom("vforce",&vforce,DOUBLE,3); + add_peratom("ervelforce",&ervelforce,DOUBLE,0); + add_peratom("etag",&etag,INT,0); + // USER-DPD package - add_peratom("rho",&eradius,DOUBLE,0); add_peratom("dpdTheta",&dpdTheta,DOUBLE,0); add_peratom("uCond",&uCond,DOUBLE,0); add_peratom("uMech",&uMech,DOUBLE,0); @@ -548,7 +567,26 @@ void Atom::peratom_create() add_peratom("edpd_cv",&edpd_cv,DOUBLE,0); add_peratom("edpd_temp",&edpd_temp,DOUBLE,0); add_peratom("edpd_flux",&edpd_flux,DOUBLE,0,1); // set per-thread flag - add_peratom("vest",&vest,DOUBLE,4); + add_peratom("cc",&cc,DOUBLE,1); + add_peratom("cc_flux",&cc_flux,DOUBLE,1,1); // set per-thread flag + + // USER-SPH package + + add_peratom("rho",&rho,DOUBLE,0); + add_peratom("drho",&drho,DOUBLE,0,1); // set per-thread flag + add_peratom("e",&e,DOUBLE,0); + add_peratom("de",&de,DOUBLE,0,1); // set per-thread flag + add_peratom("vest",&vest,DOUBLE,3); + add_peratom("cv",&cv,DOUBLE,0); + + // USER-SMD package + + add_peratom("contact_radius",&contact_radius,DOUBLE,0); + add_peratom("smd_data_9",&smd_data_9,DOUBLE,1); + add_peratom("smd_stress",&smd_stress,DOUBLE,1); + add_peratom("eff_plastic_strain",&eff_plastic_strain,DOUBLE,0); + add_peratom("eff_plastic_strain_rate",&eff_plastic_strain_rate,DOUBLE,0); + add_peratom("damage",&damage,DOUBLE,0); } /* ---------------------------------------------------------------------- @@ -579,6 +617,18 @@ void Atom::add_peratom(const char *name, void *address, nperatom++; } +/* ---------------------------------------------------------------------- + change the column count fof an existing peratom array entry + allows atom_style to specify column count as an argument + see atom_style tdpd as an example +------------------------------------------------------------------------- */ + +void Atom::add_peratom_change_columns(const char *name, int cols) +{ + for (int i = 0; i < nperatom; i++) + if (strcmp(name,peratom[i].name) == 0) peratom[i].cols = cols; +} + /* ---------------------------------------------------------------------- add info for a single per-atom array to PerAtom data struct cols = address of int variable with max columns per atom diff --git a/src/atom.h b/src/atom.h index a105b3e5b1..007d364a7b 100644 --- a/src/atom.h +++ b/src/atom.h @@ -60,6 +60,8 @@ class Atom : protected Pointers { imageint *image; double **x,**v,**f; + // charged and dipolar particles + double *rmass; double *q,**mu; @@ -69,7 +71,7 @@ class Atom : protected Pointers { double **omega,**angmom,**torque; int *ellipsoid,*line,*tri,*body; - // MOLECULE package + // molecular systems tagint *molecule; int *molindex,*molatom; @@ -101,15 +103,14 @@ class Atom : protected Pointers { // SPIN package - double **sp; - double **fm; - double **fm_long; + double **sp,**fm,**fm_long; - // USER-AWPMD and USER_EFF packages + // USER_EFF and USER-AWPMD packages int *spin; - double *eradius,*ervel,*erforce,*ervelforce; - double *cs,*csforce,*vforce; + double *eradius,*ervel,*erforce; + double *ervelforce,*cs,*csforce; + double **vforce; int *etag; // USER-DPD package @@ -261,6 +262,7 @@ class Atom : protected Pointers { void settings(class Atom *); void peratom_create(); void add_peratom(const char *, void *, int, int, int threadflag=0); + void add_peratom_change_columns(const char *, int); void add_peratom_vary(const char *, void *, int, int *, void *, int collength=0); void create_avec(const char *, int, char **, int); diff --git a/src/atom_vec_hybrid.h b/src/atom_vec_hybrid.h index 83d29112dd..7d838b7a7f 100644 --- a/src/atom_vec_hybrid.h +++ b/src/atom_vec_hybrid.h @@ -57,11 +57,6 @@ class AtomVecHybrid : public AtomVec { void pack_data_pre(int); void pack_data_post(int); - //void create_atom_post(int); - //void data_atom_post(int); - //void pack_data_pre(int); - //void pack_data_post(int); - int property_atom(char *); void pack_property_atom(int, double *, int, int); diff --git a/src/atom_vec_sphere.cpp b/src/atom_vec_sphere.cpp index 6796f9c8f2..81f1b02fdb 100644 --- a/src/atom_vec_sphere.cpp +++ b/src/atom_vec_sphere.cpp @@ -28,6 +28,7 @@ using namespace MathConst; AtomVecSphere::AtomVecSphere(LAMMPS *lmp) : AtomVec(lmp) { + mass_type = 0; molecular = 0; atom->sphere_flag = 1; @@ -51,8 +52,6 @@ AtomVecSphere::AtomVecSphere(LAMMPS *lmp) : AtomVec(lmp) fields_create = (char *) "radius rmass omega"; fields_data_atom = (char *) "id type radius rmass x"; fields_data_vel = (char *) "id v omega"; - - setup_fields(); } /* ---------------------------------------------------------------------- @@ -74,6 +73,10 @@ void AtomVecSphere::process_args(int narg, char **arg) fields_comm = (char *) "radius rmass"; fields_comm_vel = (char *) "radius rmass omega"; + + // delay setting up of fields until now + + setup_fields(); } /* ---------------------------------------------------------------------- */