From 4aecc6eaf0f901ff648a1b6561c0019fe0c1a892 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 26 Jul 2013 16:53:49 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@10430 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/PERI/fix_peri_neigh.cpp | 104 ++++-- src/PERI/fix_peri_neigh.h | 10 +- src/PERI/pair_peri_lps.cpp | 97 +---- src/PERI/pair_peri_lps.h | 1 - src/PERI/pair_peri_pmb.h | 8 +- src/PERI/pair_peri_ves.cpp | 724 ++++++++++++++++++++++++++++++++++++ src/PERI/pair_peri_ves.h | 109 ++++++ 7 files changed, 928 insertions(+), 125 deletions(-) create mode 100644 src/PERI/pair_peri_ves.cpp create mode 100644 src/PERI/pair_peri_ves.h diff --git a/src/PERI/fix_peri_neigh.cpp b/src/PERI/fix_peri_neigh.cpp index 659f47c711..a6ac8e2881 100644 --- a/src/PERI/fix_peri_neigh.cpp +++ b/src/PERI/fix_peri_neigh.cpp @@ -12,13 +12,14 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing author: Mike Parks (SNL) + Contributing authors: Mike Parks (SNL), Ezwanur Rahman, J.T. Foster (UTSA) ------------------------------------------------------------------------- */ #include "math.h" #include "fix_peri_neigh.h" #include "pair_peri_pmb.h" #include "pair_peri_lps.h" +#include "pair_peri_ves.h" #include "atom.h" #include "domain.h" #include "force.h" @@ -40,6 +41,11 @@ using namespace FixConst; FixPeriNeigh::FixPeriNeigh(LAMMPS *lmp,int narg, char **arg) : Fix(lmp, narg, arg) { + isPMB = isLPS = isVES = 0; + if (force->pair_match("peri/pmb",1)) isPMB = 1; + if (force->pair_match("peri/lps",1)) isLPS = 1; + if (force->pair_match("peri/ves",1)) isVES = 1; + restart_global = 1; restart_peratom = 1; first = 1; @@ -51,6 +57,8 @@ FixPeriNeigh::FixPeriNeigh(LAMMPS *lmp,int narg, char **arg) : maxpartner = 1; npartner = NULL; partner = NULL; + deviatorextention = NULL; + deviatorBackextention = NULL; r0 = NULL; vinter = NULL; wvolume = NULL; @@ -82,6 +90,8 @@ FixPeriNeigh::~FixPeriNeigh() memory->destroy(npartner); memory->destroy(partner); + memory->destroy(deviatorextention); + memory->destroy(deviatorBackextention); memory->destroy(r0); memory->destroy(vinter); memory->destroy(wvolume); @@ -195,12 +205,16 @@ void FixPeriNeigh::setup(int vflag) // realloc arrays with correct value for maxpartner memory->destroy(partner); + memory->destroy(deviatorextention); + memory->destroy(deviatorBackextention); memory->destroy(r0); memory->destroy(npartner); npartner = NULL; partner = NULL; - r0 = NULL; + deviatorextention = NULL; + deviatorBackextention = NULL; + r0 = NULL; grow_arrays(atom->nmax); // create partner list and r0 values from neighbor list @@ -233,7 +247,10 @@ void FixPeriNeigh::setup(int vflag) if (rsq <= cutsq[itype][jtype]) { partner[i][npartner[i]] = tag[j]; - r0[i][npartner[i]] = sqrt(rsq); + if (isVES) + deviatorextention[i][npartner[i]] = + deviatorBackextention[i][npartner[i]] = 0.0; + r0[i][npartner[i]] = sqrt(rsq); npartner[i]++; vinter[i] += vfrac[j]; } @@ -263,6 +280,8 @@ void FixPeriNeigh::setup(int vflag) double vfrac_scale; PairPeriLPS *pairlps = static_cast(anypair); PairPeriPMB *pairpmb = static_cast(anypair); + PairPeriVES *pairves = static_cast(anypair); + for (i = 0; i < nlocal; i++) { double xtmp0 = x0[i][0]; double ytmp0 = x0[i][1]; @@ -301,12 +320,16 @@ void FixPeriNeigh::setup(int vflag) (1.0 + ((delta - half_lc)/(2*half_lc) ) ); else vfrac_scale = 1.0; - if (pairpmb != NULL) // define influence function to be 1.0 - wvolume[i] += 1.0 * rsq0 * vfrac[j] * vfrac_scale; - else if (pairlps != NULL) // call the PairPeriLPS influence function + // for PMB, influence = 1.0, otherwise invoke influence function + + if (isPMB) + wvolume[i] += 1.0 * rsq0 * vfrac[j] * vfrac_scale; + else if (isLPS) wvolume[i] += pairlps->influence_function(delx0,dely0,delz0) * rsq0 * vfrac[j] * vfrac_scale; - + else if (isVES) + wvolume[i] += pairves->influence_function(delx0,dely0,delz0) * + rsq0 * vfrac[j] * vfrac_scale; } } @@ -340,14 +363,18 @@ void FixPeriNeigh::setup(int vflag) ------------------------------------------------------------------------- */ double FixPeriNeigh::memory_usage() -{ +{ int nmax = atom->nmax; int bytes = nmax * sizeof(int); bytes += nmax*maxpartner * sizeof(int); bytes += nmax*maxpartner * sizeof(double); + if (isVES) { + bytes += nmax*maxpartner * sizeof(double); + bytes += nmax*maxpartner * sizeof(double); + } bytes += nmax * sizeof(double); bytes += nmax * sizeof(double); - return bytes; + return bytes; } /* ---------------------------------------------------------------------- @@ -356,22 +383,32 @@ double FixPeriNeigh::memory_usage() void FixPeriNeigh::grow_arrays(int nmax) { - memory->grow(npartner,nmax,"peri_neigh:npartner"); - memory->grow(partner,nmax,maxpartner,"peri_neigh:partner"); - memory->grow(r0,nmax,maxpartner,"peri_neigh:r0"); - memory->grow(vinter,nmax,"peri_neigh:vinter"); - memory->grow(wvolume,nmax,"peri_neigh:wvolume"); + memory->grow(npartner,nmax,"peri_neigh:npartner"); + memory->grow(partner,nmax,maxpartner,"peri_neigh:partner"); + if (isVES) { + memory->grow(deviatorextention,nmax,maxpartner, + "peri_neigh:deviatorextention"); + memory->grow(deviatorBackextention,nmax,maxpartner, + "peri_neigh:deviatorBackextention"); + } + memory->grow(r0,nmax,maxpartner,"peri_neigh:r0"); + memory->grow(vinter,nmax,"peri_neigh:vinter"); + memory->grow(wvolume,nmax,"peri_neigh:wvolume"); } /* ---------------------------------------------------------------------- copy values within local atom-based arrays ------------------------------------------------------------------------- */ -void FixPeriNeigh::copy_arrays(int i, int j, int delflag) +void FixPeriNeigh::copy_arrays(int i, int j) { npartner[j] = npartner[i]; for (int m = 0; m < npartner[j]; m++) { partner[j][m] = partner[i][m]; + if (isVES) { + deviatorextention[j][m] = deviatorextention[i][m]; + deviatorBackextention[j][m] = deviatorBackextention[i][m]; + } r0[j][m] = r0[i][m]; } vinter[j] = vinter[i]; @@ -391,10 +428,14 @@ int FixPeriNeigh::pack_exchange(int i, double *buf) for (int n = 0; n < npartner[i]; n++) { if (partner[i][n] == 0) continue; buf[m++] = partner[i][n]; + if (isVES) { + buf[m++] = deviatorextention[i][n]; + buf[m++] = deviatorBackextention[i][n]; + } buf[m++] = r0[i][n]; } - - buf[0] = m/2; + if (isVES) buf[0] = m/4; + else buf[0] = m/2; buf[m++] = vinter[i]; buf[m++] = wvolume[i]; return m; @@ -410,7 +451,11 @@ int FixPeriNeigh::unpack_exchange(int nlocal, double *buf) npartner[nlocal] = static_cast (buf[m++]); for (int n = 0; n < npartner[nlocal]; n++) { partner[nlocal][n] = static_cast (buf[m++]); - r0[nlocal][n] = buf[m++]; + if (isVES) { + deviatorextention[nlocal][n] = buf[m++]; + deviatorBackextention[nlocal][n] = buf[m++]; + } + r0[nlocal][n] = buf[m++]; } vinter[nlocal] = buf[m++]; wvolume[nlocal] = buf[m++]; @@ -486,15 +531,21 @@ void FixPeriNeigh::restart(char *buf) int FixPeriNeigh::pack_restart(int i, double *buf) { int m = 0; - buf[m++] = 2*npartner[i] + 4; + if (isVES) buf[m++] = 4*npartner[i] + 4; + else buf[m++] = 2*npartner[i] + 4; + buf[m++] = npartner[i]; for (int n = 0; n < npartner[i]; n++) { buf[m++] = partner[i][n]; + if (isVES) { + buf[m++] = deviatorextention[i][n]; + buf[m++] = deviatorBackextention[i][n]; + } buf[m++] = r0[i][n]; } buf[m++] = vinter[i]; buf[m++] = wvolume[i]; - return m; + return m; } /* ---------------------------------------------------------------------- @@ -503,6 +554,7 @@ int FixPeriNeigh::pack_restart(int i, double *buf) void FixPeriNeigh::unpack_restart(int nlocal, int nth) { + double **extra = atom->extra; // skip to Nth set of extra values @@ -514,10 +566,14 @@ void FixPeriNeigh::unpack_restart(int nlocal, int nth) npartner[nlocal] = static_cast (extra[nlocal][m++]); for (int n = 0; n < npartner[nlocal]; n++) { partner[nlocal][n] = static_cast (extra[nlocal][m++]); + if (isVES) { + deviatorextention[nlocal][n] = extra[nlocal][m++]; + deviatorBackextention[nlocal][n] = extra[nlocal][m++]; + } r0[nlocal][n] = extra[nlocal][m++]; } vinter[nlocal] = extra[nlocal][m++]; - wvolume[nlocal] = extra[nlocal][m++]; + wvolume[nlocal] = extra[nlocal][m++]; } /* ---------------------------------------------------------------------- @@ -526,7 +582,8 @@ void FixPeriNeigh::unpack_restart(int nlocal, int nth) int FixPeriNeigh::maxsize_restart() { - return 2*maxpartner + 4; + if (isVES) return 4*maxpartner + 4; + return 2*maxpartner + 4; } /* ---------------------------------------------------------------------- @@ -535,5 +592,6 @@ int FixPeriNeigh::maxsize_restart() int FixPeriNeigh::size_restart(int nlocal) { - return 2*npartner[nlocal] + 4; + if (isVES) return 4*npartner[nlocal] + 4; + return 2*npartner[nlocal] + 4; } diff --git a/src/PERI/fix_peri_neigh.h b/src/PERI/fix_peri_neigh.h index c38aada74f..c19fb20b3a 100644 --- a/src/PERI/fix_peri_neigh.h +++ b/src/PERI/fix_peri_neigh.h @@ -28,6 +28,7 @@ class FixPeriNeigh : public Fix { friend class PairPeriPMB; friend class PairPeriPMBOMP; friend class PairPeriLPS; + friend class PairPeriVES; //NEW friend class PairPeriLPSOMP; friend class ComputeDamageAtom; @@ -42,7 +43,7 @@ class FixPeriNeigh : public Fix { double memory_usage(); void grow_arrays(int); - void copy_arrays(int, int, int); + void copy_arrays(int, int); int pack_exchange(int, double *); int unpack_exchange(int, double *); void write_restart(FILE *); @@ -60,9 +61,16 @@ class FixPeriNeigh : public Fix { int maxpartner; // max # of peridynamic neighs for any atom int *npartner; // # of neighbors for each atom int **partner; // neighs for each atom, stored as global IDs + double **deviatorextention; // Deviatoric extention + double **deviatorBackextention; // Deviatoric back extention double **r0; // initial distance to partners + double **r1; // Instanteneous distance to partners *** NEW *** + double *thetaOld; // Dilatation Old one double *vinter; // sum of vfrac for bonded neighbors double *wvolume; // weighted volume of particle + int isPMB; + int isLPS; + int isVES; class NeighList *list; }; diff --git a/src/PERI/pair_peri_lps.cpp b/src/PERI/pair_peri_lps.cpp index 6ef47b36f6..a140c77e5c 100644 --- a/src/PERI/pair_peri_lps.cpp +++ b/src/PERI/pair_peri_lps.cpp @@ -42,6 +42,7 @@ PairPeriLPS::PairPeriLPS(LAMMPS *lmp) : Pair(lmp) { for (int i = 0; i < 6; i++) virial[i] = 0.0; no_virial_fdotr_compute = 1; + single_enable = 0; ifix_peri = -1; @@ -503,102 +504,6 @@ void PairPeriLPS::read_restart(FILE *fp) } } -/* ---------------------------------------------------------------------- */ - -double PairPeriLPS::single(int i, int j, int itype, int jtype, - double rsq, double factor_coul, double factor_lj, - double &fforce) -{ - double delx0,dely0,delz0,rsq0; - double d_ij,r,dr,rk,vfrac_scale; - - double *vfrac = atom->vfrac; - double **x0 = atom->x0; - double **r0 = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0; - int **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner; - int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner; - double *wvolume = ((FixPeriNeigh *) modify->fix[ifix_peri])->wvolume; - - double lc = domain->lattice->xlattice; - double half_lc = 0.5*lc; - - double kshort; - - delx0 = x0[i][0] - x0[j][0]; - dely0 = x0[i][1] - x0[j][1]; - delz0 = x0[i][2] - x0[j][2]; - int periodic = domain->xperiodic || domain->yperiodic || domain->zperiodic; - if (periodic) domain->minimum_image(delx0,dely0,delz0); - rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0; - - d_ij = MIN(0.9*sqrt(rsq0),1.35*lc); - r = sqrt(rsq); - - double energy = 0.0; - fforce = 0.0; - - if (r < d_ij) { - dr = r - d_ij; - // kshort resembles short-range force constant of bond-based theory in 3d - kshort = (15.0 * 18.0 * bulkmodulus[itype][itype]) / - ( 3.141592653589793 * cutsq[itype][jtype] * cutsq[itype][jtype]); - rk = ( kshort * vfrac[j]) * (dr / sqrt(cutsq[itype][jtype])); - if (r > 0.0) fforce += -(rk/r); - energy += 0.5*rk*dr; - } - - if (atom->nmax > nmax) { - memory->destroy(theta); - nmax = atom->nmax; - memory->create(theta,nmax,"pair:theta"); - } - - // Compute the dilatation on each particle - compute_dilatation(); - - // communicate dilatation (theta) of each particle - comm->forward_comm_pair(this); - // communicate wighted volume (wvolume) upon every reneighbor - if (neighbor->ago == 0) - comm->forward_comm_fix(modify->fix[ifix_peri]); - - double omega_plus, omega_minus; - - int jnum = npartner[i]; - for (int jj = 0; jj < jnum; jj++) { - if (partner[i][jj] == 0) continue; - if (j < 0) continue; - if (j == atom->map(partner[i][jj])) { - dr = r - r0[i][jj]; - if (fabs(dr) < 2.2204e-016) dr = 0.0; - - // scale vfrac[j] if particle j near the horizon - - if ( (fabs(r0[i][jj] - sqrt(cutsq[itype][jtype]))) <= half_lc) - vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) + - (1.0 + ((sqrt(cutsq[itype][jtype]) - half_lc)/(2*half_lc))); - else vfrac_scale = 1.0; - - omega_plus = influence_function(-1.0*delx0,-1.0*dely0,-1.0*delz0); - omega_minus = influence_function(delx0,dely0,delz0); - rk = (3.0* bulkmodulus[itype][itype] -5.0 * shearmodulus[itype][itype]) * - vfrac[j] * vfrac_scale * ( (omega_plus * theta[i] / wvolume[i]) + - (omega_minus * theta[j] / wvolume[j])) * - r0[i][jj]; - rk += 15.0 * ( shearmodulus[itype][itype] * vfrac[j] * vfrac_scale ) * - ( (omega_plus / wvolume[i]) + (omega_minus / wvolume[j]) ) * dr; - - if (r > 0.0) fforce += -(rk/r); - energy += 0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) * - omega_plus * ( dr - theta[i]* r0[i][jj] / 3.0 ) * - ( dr - theta[i]* r0[i][jj] / 3.0 ) * vfrac[j] * vfrac_scale; - - } - } - - return energy; -} - /* ---------------------------------------------------------------------- memory usage of local atom-based arrays ------------------------------------------------------------------------- */ diff --git a/src/PERI/pair_peri_lps.h b/src/PERI/pair_peri_lps.h index 22d901e4f8..cadbca6dd2 100644 --- a/src/PERI/pair_peri_lps.h +++ b/src/PERI/pair_peri_lps.h @@ -40,7 +40,6 @@ class PairPeriLPS : public Pair { void read_restart(FILE *); void write_restart_settings(FILE *) {} void read_restart_settings(FILE *) {} - double single(int, int, int, int, double, double, double, double &); double memory_usage(); double influence_function(double, double, double); void compute_dilatation(); diff --git a/src/PERI/pair_peri_pmb.h b/src/PERI/pair_peri_pmb.h index b69fe4c3a7..81e78bb368 100644 --- a/src/PERI/pair_peri_pmb.h +++ b/src/PERI/pair_peri_pmb.h @@ -83,6 +83,10 @@ E: Pair peri requires an atom map, see atom_modify Even for atomic systems, an atom map is required to find Peridynamic bonds. Use the atom_modify command to define one. +E: Pair peri requires a lattice be defined + +Use the lattice command for this purpose. + E: Pair peri lattice is not identical in x, y, and z The lattice defined by the lattice command must be cubic. @@ -91,8 +95,4 @@ E: Fix peri neigh does not exist Somehow a fix that the pair style defines has been deleted. -U: Pair peri requires a lattice be defined - -Use the lattice command for this purpose. - */ diff --git a/src/PERI/pair_peri_ves.cpp b/src/PERI/pair_peri_ves.cpp new file mode 100644 index 0000000000..5586abc965 --- /dev/null +++ b/src/PERI/pair_peri_ves.cpp @@ -0,0 +1,724 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing authors: Ezwanur Rahman, J.T. Foster (UTSA) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdlib.h" +#include "string.h" +#include "pair_peri_ves.h" +#include "atom.h" +#include "domain.h" +#include "lattice.h" +#include "force.h" +#include "update.h" +#include "modify.h" +#include "fix.h" +#include "fix_peri_neigh.h" +#include "comm.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "memory.h" +#include "error.h" +#include "update.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +PairPeriVES::PairPeriVES(LAMMPS *lmp) : Pair(lmp) +{ + for (int i = 0; i < 6; i++) virial[i] = 0.0; + no_virial_fdotr_compute = 1; + single_enable = 0; + + ifix_peri = -1; + + nmax = 0; + s0_new = NULL; + theta = NULL; + + bulkmodulus = NULL; + shearmodulus = NULL; + s00 = alpha = NULL; + cut = NULL; + m_lambdai = NULL; + m_taubi = NULL; + + // set comm size needed by this Pair + // comm_reverse not needed + + comm_forward = 1; +} + +/* ---------------------------------------------------------------------- */ + +PairPeriVES::~PairPeriVES() +{ + if (ifix_peri >= 0) modify->delete_fix("PERI_NEIGH"); + + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + memory->destroy(bulkmodulus); + memory->destroy(shearmodulus); + memory->destroy(s00); + memory->destroy(alpha); + memory->destroy(cut); + memory->destroy(m_lambdai); + memory->destroy(m_taubi); + memory->destroy(theta); + memory->destroy(s0_new); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairPeriVES::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz; + double xtmp0,ytmp0,ztmp0,delx0,dely0,delz0,rsq0; + double rsq,r,dr,dr1,rk,evdwl,fpair,fbond; + double deltaed,fbondViscoElastic,fbondFinal; + double decay,betai,lambdai,edbNp1,rkNew; + int *ilist,*jlist,*numneigh,**firstneigh; + double d_ij,delta,stretch; + + evdwl = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = eflag_global = eflag_atom = 0; + + double **f = atom->f; + double **x = atom->x; + int *type = atom->type; + int nlocal = atom->nlocal; + + double timestepsize = update->dt; + double *vfrac = atom->vfrac; + double *s0 = atom->s0; + double **x0 = atom->x0; + double **r0 = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0; + double **deviatorextention = + ((FixPeriNeigh *) modify->fix[ifix_peri])->deviatorextention; + double **deviatorBackextention = + ((FixPeriNeigh *) modify->fix[ifix_peri])->deviatorBackextention; + int **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner; + int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner; + double *wvolume = ((FixPeriNeigh *) modify->fix[ifix_peri])->wvolume; + + // lc = lattice constant + // init_style guarantees it's the same in x, y, and z + + double lc = domain->lattice->xlattice; + double half_lc = 0.5*lc; + double vfrac_scale = 1.0; + + // short-range forces + + int newton_pair = force->newton_pair; + int periodic = domain->xperiodic || domain->yperiodic || domain->zperiodic; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + // need minimg() for x0 difference since not ghosted + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + xtmp0 = x0[i][0]; + ytmp0 = x0[i][1]; + ztmp0 = x0[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + + rsq = delx*delx + dely*dely + delz*delz; + delx0 = xtmp0 - x0[j][0]; + dely0 = ytmp0 - x0[j][1]; + delz0 = ztmp0 - x0[j][2]; + if (periodic) domain->minimum_image(delx0,dely0,delz0); + rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0; + jtype = type[j]; + + r = sqrt(rsq); + + // short-range interaction distance based on initial particle position + // 0.9 and 1.35 are constants + + d_ij = MIN(0.9*sqrt(rsq0),1.35*lc); + + // short-range contact forces + // 15 is constant taken from the EMU Theory Manual + // Silling, 12 May 2005, p 18 + + if (r < d_ij) { + dr = r - d_ij; + + // kshort based upon short-range force constant + // of the bond-based theory used in PMB model + + double kshort = (15.0 * 18.0 * bulkmodulus[itype][itype]) / + (3.141592653589793 * cutsq[itype][jtype] * cutsq[itype][jtype]); + rk = (kshort * vfrac[j]) * (dr / cut[itype][jtype]); + + if (r > 0.0) fpair = -(rk/r); + else fpair = 0.0; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if (eflag) evdwl = 0.5*rk*dr; + if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0, + fpair*vfrac[i],delx,dely,delz); + } + } + } + + // grow bond forces array if necessary + + if (atom->nmax > nmax) { + memory->destroy(s0_new); + memory->destroy(theta); + nmax = atom->nmax; + memory->create(s0_new,nmax,"pair:s0_new"); + memory->create(theta,nmax,"pair:theta"); + } + + // compute the dilatation on each particle + + compute_dilatation(); + + // communicate dilatation (theta) of each particle + + comm->forward_comm_pair(this); + + // communicate weighted volume (wvolume) upon every reneighbor + + if (neighbor->ago == 0) + comm->forward_comm_fix(modify->fix[ifix_peri]); + + // volume-dependent part of the energy + + if (eflag) { + for (i = 0; i < nlocal; i++) { + itype = type[i]; + if (eflag_global) + eng_vdwl += 0.5 * bulkmodulus[itype][itype] * (theta[i] * theta[i]); + if (eflag_atom) + eatom[i] += 0.5 * bulkmodulus[itype][itype] * (theta[i] * theta[i]); + } + } + + // loop over my particles and their partners + // partner list contains all bond partners, so I-J appears twice + // if bond already broken, skip this partner + // first = true if this is first neighbor of particle i + + bool first; + double omega_minus, omega_plus; + + for (i = 0; i < nlocal; i++) { + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + xtmp0 = x0[i][0]; + ytmp0 = x0[i][1]; + ztmp0 = x0[i][2]; + itype = type[i]; + jnum = npartner[i]; + first = true; + + for (jj = 0; jj < jnum; jj++) { + if (partner[i][jj] == 0) continue; + j = atom->map(partner[i][jj]); + + // check if lost a partner without first breaking bond + + if (j < 0) { + partner[i][jj] = 0; + continue; + } + + // compute force density, add to PD equation of motion + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + if (periodic) domain->minimum_image(delx,dely,delz); + rsq = delx*delx + dely*dely + delz*delz; + delx0 = xtmp0 - x0[j][0]; + dely0 = ytmp0 - x0[j][1]; + delz0 = ztmp0 - x0[j][2]; + if (periodic) domain->minimum_image(delx0,dely0,delz0); + jtype = type[j]; + delta = cut[itype][jtype]; + r = sqrt(rsq); + dr = r - r0[i][jj]; + + // avoid roundoff errors + + if (fabs(dr) < 2.2204e-016) dr = 0.0; + + // scale vfrac[j] if particle j near the horizon + + if ((fabs(r0[i][jj] - delta)) <= half_lc) + vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) + + (1.0 + ((delta - half_lc)/(2*half_lc) ) ); + else vfrac_scale = 1.0; + + omega_plus = influence_function(-1.0*delx0,-1.0*dely0,-1.0*delz0); + omega_minus = influence_function(delx0,dely0,delz0); + + rk = ( (3.0 * bulkmodulus[itype][itype]) * vfrac[j] * vfrac_scale * + ( (omega_plus * theta[i] / wvolume[i]) + + ( omega_minus * theta[j] / wvolume[j] ) ) ) * r0[i][jj]; + + if (r > 0.0) fbond = -(rk/r); + else fbond = 0.0; + + // for viscoelasticity + + lambdai=m_lambdai[itype][itype]; + double taui = m_taubi[itype][itype]; + double c1 = taui/timestepsize; + decay=exp(-1.0/c1); + betai=1.-c1*(1.-decay); + + double deviatoric_extension = + dr - (theta[i]* r0[i][jj] / 3.0); + deltaed = deviatoric_extension-deviatorextention[i][jj]; + + // back extention at current step + + edbNp1 = deviatorextention[i][jj]*(1-decay) + + deviatorBackextention[i][jj]*decay+betai*deltaed; + + rkNew = ((1-lambdai)*15.0) * + ( shearmodulus[itype][itype] * vfrac[j] * vfrac_scale ) * + ( (omega_plus / wvolume[i]) + (omega_minus / wvolume[j]) ) * + deviatoric_extension; + rkNew += (lambdai*15.0) * + ( shearmodulus[itype][itype] * vfrac[j] * vfrac_scale ) * + ( (omega_plus / wvolume[i]) + (omega_minus / wvolume[j]) ) * + (deviatoric_extension-edbNp1); + + if (r > 0.0) fbondViscoElastic = -(rkNew/r); + else fbondViscoElastic = 0.0; + + // total Force: elastic + viscoelastic + + fbondFinal=fbond+fbondViscoElastic; + fbond=fbondFinal; + + f[i][0] += delx*fbond; + f[i][1] += dely*fbond; + f[i][2] += delz*fbond; + + // since I-J is double counted, set newton off & use 1/2 factor and I,I + + if (eflag) evdwl = (0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) * + omega_plus * deviatoric_extension * + deviatoric_extension) + + (0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) * + omega_plus * (deviatoric_extension-edbNp1) * + (deviatoric_extension-edbNp1)) * vfrac[j] * vfrac_scale; + if (evflag) ev_tally(i,i,nlocal,0,0.5*evdwl,0.0, + 0.5*fbond*vfrac[i],delx,dely,delz); + + // find stretch in bond I-J and break if necessary + // use s0 from previous timestep + + // store current deviatoric extention + + deviatorextention[i][jj]=deviatoric_extension; + deviatorBackextention[i][jj]=edbNp1; + + stretch = dr / r0[i][jj]; + if (stretch > MIN(s0[i],s0[j])) partner[i][jj] = 0; + + // update s0 for next timestep + + if (first) + s0_new[i] = s00[itype][jtype] - (alpha[itype][jtype] * stretch); + else + s0_new[i] = MAX(s0_new[i],s00[itype][jtype] - + (alpha[itype][jtype] * stretch)); + + first = false; + } + } + + // store new s0 + + for (i = 0; i < nlocal; i++) s0[i] = s0_new[i]; +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairPeriVES::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + memory->create(setflag,n+1,n+1,"pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + memory->create(cutsq,n+1,n+1,"pair:cutsq"); + memory->create(bulkmodulus,n+1,n+1,"pair:bulkmodulus"); + memory->create(shearmodulus,n+1,n+1,"pair:shearmodulus"); + memory->create(s00,n+1,n+1,"pair:s00"); + memory->create(alpha,n+1,n+1,"pair:alpha"); + memory->create(cut,n+1,n+1,"pair:cut"); + memory->create(m_lambdai,n+1,n+1,"pair:m_lambdai"); + memory->create(m_taubi,n+1,n+1,"pair:m_taubi"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairPeriVES::settings(int narg, char **arg) +{ + if (narg) error->all(FLERR,"Illegal pair_style command"); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairPeriVES::coeff(int narg, char **arg) +{ + if (narg != 9) error->all(FLERR,"Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + force->bounds(arg[0],atom->ntypes,ilo,ihi); + force->bounds(arg[1],atom->ntypes,jlo,jhi); + + double bulkmodulus_one = atof(arg[2]); + double shearmodulus_one = atof(arg[3]); + double cut_one = atof(arg[4]); + double s00_one = atof(arg[5]); + double alpha_one = atof(arg[6]); + double mlambdai_one = atof(arg[7]); + double mtaui_one = atof(arg[8]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + bulkmodulus[i][j] = bulkmodulus_one; + shearmodulus[i][j] = shearmodulus_one; + cut[i][j] = cut_one; + s00[i][j] = s00_one; + alpha[i][j] = alpha_one; + m_lambdai[i][j] = mlambdai_one; + m_taubi[i][j] = mtaui_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairPeriVES::init_one(int i, int j) +{ + if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); + + bulkmodulus[j][i] = bulkmodulus[i][j]; + shearmodulus[j][i] = shearmodulus[i][j]; + s00[j][i] = s00[i][j]; + alpha[j][i] = alpha[i][j]; + cut[j][i] = cut[i][j]; + m_lambdai[j][i] = m_lambdai[i][j]; + m_taubi[j][i] = m_taubi[i][j]; + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairPeriVES::init_style() +{ + // error checks + + if (!atom->peri_flag) + error->all(FLERR,"Pair style peri requires atom style peri"); + if (atom->map_style == 0) + error->all(FLERR,"Pair peri requires an atom map, see atom_modify"); + + if (domain->lattice == NULL) + error->all(FLERR,"Pair peri requires a lattice be defined"); + if (domain->lattice->xlattice != domain->lattice->ylattice || + domain->lattice->xlattice != domain->lattice->zlattice || + domain->lattice->ylattice != domain->lattice->zlattice) + error->all(FLERR,"Pair peri lattice is not identical in x, y, and z"); + + // if first init, create Fix needed for storing fixed neighbors + + if (ifix_peri == -1) { + char **fixarg = new char*[3]; + fixarg[0] = (char *) "PERI_NEIGH"; + fixarg[1] = (char *) "all"; + fixarg[2] = (char *) "PERI_NEIGH"; + modify->add_fix(3,fixarg); + delete [] fixarg; + } + + // find associated PERI_NEIGH fix that must exist + // could have changed locations in fix list since created + + for (int i = 0; i < modify->nfix; i++) + if (strcmp(modify->fix[i]->style,"PERI_NEIGH") == 0) ifix_peri = i; + if (ifix_peri == -1) error->all(FLERR,"Fix peri neigh does not exist"); + + neighbor->request(this); +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairPeriVES::write_restart(FILE *fp) +{ + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + fwrite(&setflag[i][j],sizeof(int),1,fp); + if (setflag[i][j]) { + fwrite(&bulkmodulus[i][j],sizeof(double),1,fp); + fwrite(&shearmodulus[i][j],sizeof(double),1,fp); + fwrite(&s00[i][j],sizeof(double),1,fp); + fwrite(&alpha[i][j],sizeof(double),1,fp); + fwrite(&cut[i][j],sizeof(double),1,fp); + fwrite(&m_lambdai[i][j],sizeof(double),1,fp); + fwrite(&m_taubi[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairPeriVES::read_restart(FILE *fp) +{ + allocate(); + + int i,j; + int me = comm->me; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); + MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (setflag[i][j]) { + if (me == 0) { + fread(&bulkmodulus[i][j],sizeof(double),1,fp); + fread(&shearmodulus[i][j],sizeof(double),1,fp); + fread(&s00[i][j],sizeof(double),1,fp); + fread(&alpha[i][j],sizeof(double),1,fp); + fread(&cut[i][j],sizeof(double),1,fp); + fread(&m_lambdai[i][j],sizeof(double),1,fp); + fread(&m_taubi[i][j],sizeof(double),1,fp); + } + MPI_Bcast(&bulkmodulus[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&shearmodulus[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&s00[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&alpha[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&m_lambdai[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&m_taubi[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based arrays +------------------------------------------------------------------------- */ + +double PairPeriVES::memory_usage() +{ + double bytes = 2 * nmax * sizeof(double); + return bytes; +} + +/* ---------------------------------------------------------------------- + influence function definition +------------------------------------------------------------------------- */ + +double PairPeriVES::influence_function(double xi_x, double xi_y, double xi_z) +{ + double r = sqrt(xi_x*xi_x + xi_y*xi_y + xi_z*xi_z); + double omega; + + if (fabs(r) < 2.2204e-016) + error->one(FLERR,"Divide by 0 in influence function of pair peri/lps"); + omega = 1.0/r; + return omega; +} + +/* ---------------------------------------------------------------------- */ + +void PairPeriVES::compute_dilatation() +{ + int i,j,jj,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz; + double xtmp0,ytmp0,ztmp0,delx0,dely0,delz0; + double rsq,r,dr; + double delta; + + double **x = atom->x; + int *type = atom->type; + double **x0 = atom->x0; + int nlocal = atom->nlocal; + double *vfrac = atom->vfrac; + double vfrac_scale = 1.0; + + double lc = domain->lattice->xlattice; + double half_lc = 0.5*lc; + + double **r0 = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0; + int **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner; + int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner; + double *wvolume = ((FixPeriNeigh *) modify->fix[ifix_peri])->wvolume; + + int periodic = domain->xperiodic || domain->yperiodic || domain->zperiodic; + + // compute the dilatation theta + + for (i = 0; i < nlocal; i++) { + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + xtmp0 = x0[i][0]; + ytmp0 = x0[i][1]; + ztmp0 = x0[i][2]; + jnum = npartner[i]; + theta[i] = 0.0; + itype = type[i]; + + for (jj = 0; jj < jnum; jj++) { + + // if bond already broken, skip this partner + + if (partner[i][jj] == 0) continue; + + // look up local index of this partner particle + + j = atom->map(partner[i][jj]); + + // skip if particle is "lost" + + if (j < 0) continue; + + // compute force density and add to PD equation of motion + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + if (periodic) domain->minimum_image(delx,dely,delz); + rsq = delx*delx + dely*dely + delz*delz; + delx0 = xtmp0 - x0[j][0]; + dely0 = ytmp0 - x0[j][1]; + delz0 = ztmp0 - x0[j][2]; + if (periodic) domain->minimum_image(delx0,dely0,delz0); + + r = sqrt(rsq); + dr = r - r0[i][jj]; + if (fabs(dr) < 2.2204e-016) dr = 0.0; + + jtype = type[j]; + delta = cut[itype][jtype]; + + // scale vfrac[j] if particle j near the horizon + + if ((fabs(r0[i][jj] - delta)) <= half_lc) + vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) + + (1.0 + ((delta - half_lc)/(2*half_lc) ) ); + else vfrac_scale = 1.0; + + theta[i] += influence_function(delx0, dely0, delz0) * r0[i][jj] * dr * + vfrac[j] * vfrac_scale; + + } + + // if wvolume[i] is zero, then particle i has no bonds + // therefore, the dilatation is set to + + if (wvolume[i] != 0.0) theta[i] = (3.0/wvolume[i]) * theta[i]; + else theta[i] = 0; + } +} + +/* ---------------------------------------------------------------------- + communication routines +---------------------------------------------------------------------- */ + +int PairPeriVES::pack_comm(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = theta[j]; + } + return 1; +} + +/* ---------------------------------------------------------------------- */ + +void PairPeriVES::unpack_comm(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + theta[i] = buf[m++]; + } +} diff --git a/src/PERI/pair_peri_ves.h b/src/PERI/pair_peri_ves.h new file mode 100644 index 0000000000..c21e60c9b8 --- /dev/null +++ b/src/PERI/pair_peri_ves.h @@ -0,0 +1,109 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(peri/ves,PairPeriVES) + +#else + +#ifndef LMP_PAIR_PERI_VES_H +#define LMP_PAIR_PERI_VES_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairPeriVES : public Pair { + public: + PairPeriVES(class LAMMPS *); + virtual ~PairPeriVES(); + int pack_comm(int, int *, double *, int, int *); + void unpack_comm(int, int, double *); + + virtual void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + double init_one(int, int); + void init_style(); + void write_restart(FILE *); + void read_restart(FILE *); + void write_restart_settings(FILE *) {} + void read_restart_settings(FILE *) {} + double memory_usage(); + double influence_function(double, double, double); + void compute_dilatation(); + + protected: + int ifix_peri; + double **bulkmodulus; + double **shearmodulus; + double **s00, **alpha; + double **cut, **m_lambdai, **m_taubi; //NEW: **m_lambdai,**m_taubi + + double *s0_new; + double *theta; + int nmax; + + void allocate(); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Incorrect args for pair coefficients + +Self-explanatory. Check the input script or data file. + +E: All pair coeffs are not set + +All pair coefficients must be set in the data file or by the +pair_coeff command before running a simulation. + +E: Pair style peri requires atom style peri + +Self-explanatory. + +E: Pair peri requires an atom map, see atom_modify + +Even for atomic systems, an atom map is required to find Peridynamic +bonds. Use the atom_modify command to define one. + +E: Pair peri requires a lattice be defined + +Use the lattice command for this purpose. + +E: Pair peri lattice is not identical in x, y, and z + +The lattice defined by the lattice command must be cubic. + +E: Fix peri neigh does not exist + +Somehow a fix that the pair style defines has been deleted. + +E: Divide by 0 in influence function of pair peri/lps + +This should not normally occur. It is likely a problem with your +model. + +*/