From c6ab9c0ce6142cb228110c982680f462081d3de9 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 26 Oct 2011 17:08:27 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@7209 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/FLD/Install.sh | 35 + src/FLD/pair_brownian.cpp | 611 ++++++++++ src/FLD/pair_brownian.h | 58 + src/FLD/pair_brownian_poly.cpp | 338 ++++++ src/FLD/pair_brownian_poly.h | 39 + src/FLD/pair_lubricate.cpp | 723 +++++++++++ src/FLD/pair_lubricate.h | 59 + src/FLD/pair_lubricateU.cpp | 1933 ++++++++++++++++++++++++++++++ src/FLD/pair_lubricateU.h | 74 ++ src/FLD/pair_lubricateU_poly.cpp | 1134 ++++++++++++++++++ src/FLD/pair_lubricateU_poly.h | 45 + src/FLD/pair_lubricate_poly.cpp | 459 +++++++ src/FLD/pair_lubricate_poly.h | 38 + src/Makefile | 2 +- 14 files changed, 5547 insertions(+), 1 deletion(-) create mode 100644 src/FLD/Install.sh create mode 100755 src/FLD/pair_brownian.cpp create mode 100644 src/FLD/pair_brownian.h create mode 100644 src/FLD/pair_brownian_poly.cpp create mode 100644 src/FLD/pair_brownian_poly.h create mode 100755 src/FLD/pair_lubricate.cpp create mode 100644 src/FLD/pair_lubricate.h create mode 100644 src/FLD/pair_lubricateU.cpp create mode 100644 src/FLD/pair_lubricateU.h create mode 100644 src/FLD/pair_lubricateU_poly.cpp create mode 100644 src/FLD/pair_lubricateU_poly.h create mode 100644 src/FLD/pair_lubricate_poly.cpp create mode 100644 src/FLD/pair_lubricate_poly.h diff --git a/src/FLD/Install.sh b/src/FLD/Install.sh new file mode 100644 index 0000000000..362dad34f3 --- /dev/null +++ b/src/FLD/Install.sh @@ -0,0 +1,35 @@ +# Install/unInstall package files in LAMMPS + +if (test $1 == 1) then + + cp pair_brownian.cpp .. + cp pair_brownian_poly.cpp .. + cp pair_lubricate.cpp .. + cp pair_lubricate_poly.cpp .. + cp pair_lubricateU.cpp .. + cp pair_lubricateU_poly.cpp .. + + cp pair_brownian.h .. + cp pair_brownian_poly.h .. + cp pair_lubricate.h .. + cp pair_lubricate_poly.h .. + cp pair_lubricateU.h .. + cp pair_lubricateU_poly.h .. + +elif (test $1 == 0) then + + rm ../pair_brownian.cpp + rm ../pair_brownian_poly.cpp + rm ../pair_lubricate.cpp + rm ../pair_lubricate_poly.cpp + rm ../pair_lubricateU.cpp + rm ../pair_lubricateU_poly.cpp + + rm ../pair_brownian.h + rm ../pair_brownian_poly.h + rm ../pair_lubricate.h + rm ../pair_lubricate_poly.h + rm ../pair_lubricateU.h + rm ../pair_lubricateU_poly.h + +fi diff --git a/src/FLD/pair_brownian.cpp b/src/FLD/pair_brownian.cpp new file mode 100755 index 0000000000..61882de036 --- /dev/null +++ b/src/FLD/pair_brownian.cpp @@ -0,0 +1,611 @@ +/* ---------------------------------------------------------------------- + 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: Amit Kumar and Michael Bybee (UIUC) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_brownian.h" +#include "atom.h" +#include "atom_vec.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "domain.h" +#include "update.h" +#include "memory.h" +#include "random_mars.h" +#include "math_const.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +/* ---------------------------------------------------------------------- */ + +PairBrownian::PairBrownian(LAMMPS *lmp) : Pair(lmp) +{ + single_enable = 0; + random = NULL; +} + +/* ---------------------------------------------------------------------- */ + +PairBrownian::~PairBrownian() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + memory->destroy(cut); + memory->destroy(cut_inner); + } + delete random; +} + +/* ---------------------------------------------------------------------- */ + +void PairBrownian::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,tx,ty,tz; + double rsq,r,h_sep,radi; + int *ilist,*jlist,*numneigh,**firstneigh; + + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + double **torque = atom->torque; + double *radius = atom->radius; + int *type = atom->type; + int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; + + double vxmu2f = force->vxmu2f; + double randr; + double prethermostat; + double xl[3],a_sq,a_sh,a_pu,Fbmag; + double p1[3],p2[3],p3[3]; + int overlaps = 0; + + // scale factor for Brownian moments + + prethermostat = sqrt(24.0*force->boltz*t_target/update->dt); + prethermostat *= sqrt(force->vxmu2f/force->ftm2v/force->mvv2e); + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + radi = radius[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + // FLD contribution to force and torque due to isotropic terms + + if (flagfld) { + f[i][0] += prethermostat*sqrt(R0)*(random->uniform()-0.5); + f[i][1] += prethermostat*sqrt(R0)*(random->uniform()-0.5); + f[i][2] += prethermostat*sqrt(R0)*(random->uniform()-0.5); + if (flaglog) { + torque[i][0] += prethermostat*sqrt(RT0)*(random->uniform()-0.5); + torque[i][1] += prethermostat*sqrt(RT0)*(random->uniform()-0.5); + torque[i][2] += prethermostat*sqrt(RT0)*(random->uniform()-0.5); + } + } + + 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; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + r = sqrt(rsq); + + // scalar resistances a_sq and a_sh + + h_sep = r - 2.0*radi; + + // check for overlaps + + if (h_sep < 0.0) overlaps++; + + // if less than minimum gap, use minimum gap instead + + if (r < cut_inner[itype][jtype]) + h_sep = cut_inner[itype][jtype] - 2.0*radi; + + // scale h_sep by radi + + h_sep = h_sep/radi; + + // scalar resistances + + if (flaglog) { + a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep + 9.0/40.0*log(1.0/h_sep)); + a_sh = 6.0*MY_PI*mu*radi*(1.0/6.0*log(1.0/h_sep)); + a_pu = 8.0*MY_PI*mu*pow(radi,3)*(3.0/160.0*log(1.0/h_sep)); + } else + a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep); + + // generate the Pairwise Brownian Force: a_sq + + Fbmag = prethermostat*sqrt(a_sq); + + // generate a random number + + randr = random->uniform()-0.5; + + // contribution due to Brownian motion + + fx = Fbmag*randr*delx/r; + fy = Fbmag*randr*dely/r; + fz = Fbmag*randr*delz/r; + + // add terms due to a_sh + + if (flaglog) { + + // generate two orthogonal vectors to the line of centers + + p1[0] = delx/r; p1[1] = dely/r; p1[2] = delz/r; + set_3_orthogonal_vectors(p1,p2,p3); + + // magnitude + + Fbmag = prethermostat*sqrt(a_sh); + + // force in each of the two directions + + randr = random->uniform()-0.5; + fx += Fbmag*randr*p2[0]; + fy += Fbmag*randr*p2[1]; + fz += Fbmag*randr*p2[2]; + + randr = random->uniform()-0.5; + fx += Fbmag*randr*p3[0]; + fy += Fbmag*randr*p3[1]; + fz += Fbmag*randr*p3[2]; + } + + // scale forces to appropriate units + + fx = vxmu2f*fx; + fy = vxmu2f*fy; + fz = vxmu2f*fz; + + // sum to total force + + f[i][0] -= fx; + f[i][1] -= fy; + f[i][2] -= fz; + + if (newton_pair || j < nlocal) { + //randr = random->uniform()-0.5; + //fx = Fbmag*randr*delx/r; + //fy = Fbmag*randr*dely/r; + //fz = Fbmag*randr*delz/r; + + f[j][0] += fx; + f[j][1] += fy; + f[j][2] += fz; + } + + // torque due to the Brownian Force + + if (flaglog) { + + // location of the point of closest approach on I from its center + + xl[0] = -delx/r*radi; + xl[1] = -dely/r*radi; + xl[2] = -delz/r*radi; + + // torque = xl_cross_F + + tx = xl[1]*fz - xl[2]*fy; + ty = xl[2]*fx - xl[0]*fz; + tz = xl[0]*fy - xl[1]*fx; + + // torque is same on both particles + + torque[i][0] -= tx; + torque[i][1] -= ty; + torque[i][2] -= tz; + + if (newton_pair || j < nlocal) { + torque[j][0] -= tx; + torque[j][1] -= ty; + torque[j][2] -= tz; + } + + // torque due to a_pu + + Fbmag = prethermostat*sqrt(a_pu); + + // force in each direction + + randr = random->uniform()-0.5; + tx = Fbmag*randr*p2[0]; + ty = Fbmag*randr*p2[1]; + tz = Fbmag*randr*p2[2]; + + randr = random->uniform()-0.5; + tx += Fbmag*randr*p3[0]; + ty += Fbmag*randr*p3[1]; + tz += Fbmag*randr*p3[2]; + + // torque has opposite sign on two particles + + torque[i][0] -= tx; + torque[i][1] -= ty; + torque[i][2] -= tz; + + if (newton_pair || j < nlocal) { + torque[j][0] += tx; + torque[j][1] += ty; + torque[j][2] += tz; + } + } + + if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair, + 0.0,0.0,-fx,-fy,-fz,delx,dely,delz); + } + } + } + + int print_overlaps = 0; + if (print_overlaps && overlaps) + printf("Number of overlaps=%d\n",overlaps); + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairBrownian::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(cut,n+1,n+1,"pair:cut"); + memory->create(cut_inner,n+1,n+1,"pair:cut_inner"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairBrownian::settings(int narg, char **arg) +{ + if (narg != 7) error->all(FLERR,"Illegal pair_style command"); + + mu = atof(arg[0]); + flaglog = atoi(arg[1]); + flagfld = atoi(arg[2]); + cut_inner_global = atof(arg[3]); + cut_global = atof(arg[4]); + t_target = atof(arg[5]); + seed = atoi(arg[6]); + + // initialize Marsaglia RNG with processor-unique seed + + delete random; + random = new RanMars(lmp,seed + comm->me); + + // reset cutoffs that have been explicitly set + + if (allocated) { + for (int i = 1; i <= atom->ntypes; i++) + for (int j = i+1; j <= atom->ntypes; j++) + if (setflag[i][j]) { + cut_inner[i][j] = cut_inner_global; + cut[i][j] = cut_global; + } + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairBrownian::coeff(int narg, char **arg) +{ + if (narg != 2 && narg != 4) + 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 cut_inner_one = cut_inner_global; + double cut_one = cut_global; + + if (narg == 4) { + cut_inner_one = atof(arg[2]); + cut_one = atof(arg[3]); + } + + int count = 0; + for (int i = ilo; i <= ihi; i++) + for (int j = MAX(jlo,i); j <= jhi; j++) { + cut_inner[i][j] = cut_inner_one; + cut[i][j] = cut_one; + setflag[i][j] = 1; + count++; + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairBrownian::init_style() +{ + if (!atom->sphere_flag) + error->all(FLERR,"Pair brownian requires atom style sphere"); + + // if newton off, forces between atoms ij will be double computed + // using different random numbers + + if (force->newton_pair == 0 && comm->me == 0) + error->warning(FLERR, + "Pair brownian needs newton pair on for " + "momentum conservation"); + + int irequest = neighbor->request(this); + + // insure all particles are finite-size + // for pair hybrid, should limit test to types using the pair style + + double *radius = atom->radius; + int *type = atom->type; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (radius[i] == 0.0) + error->one(FLERR,"Pair brownian requires extended particles"); + + // require monodisperse system with same radii for all types + + double rad,radtype; + for (int i = 1; i <= atom->ntypes; i++) { + if (!atom->radius_consistency(i,radtype)) + error->all(FLERR,"Pair brownian requires monodisperse particles"); + if (i > 1 && radtype != rad) + error->all(FLERR,"Pair brownian requires monodisperse particles"); + rad = radtype; + } + + // set the isotropic constants that depend on the volume fraction + // vol_T = total volume + + double vol_T = domain->xprd*domain->yprd*domain->zprd; + + // vol_P = volume of particles, assuming mono-dispersity + // vol_f = volume fraction + + double vol_P = atom->natoms*(4.0/3.0)*MY_PI*pow(rad,3); + double vol_f = vol_P/vol_T; + + // set isotropic constants + + if (flaglog == 0) { + R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f); + RT0 = 8*MY_PI*mu*pow(rad,3); // not actually needed + } else { + R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f); + RT0 = 8*MY_PI*mu*pow(rad,3)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); + } +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairBrownian::init_one(int i, int j) +{ + if (setflag[i][j] == 0) { + cut_inner[i][j] = mix_distance(cut_inner[i][i],cut_inner[j][j]); + cut[i][j] = mix_distance(cut[i][i],cut[j][j]); + } + + cut_inner[j][i] = cut_inner[i][j]; + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairBrownian::write_restart(FILE *fp) +{ + write_restart_settings(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(&cut_inner[i][j],sizeof(double),1,fp); + fwrite(&cut[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairBrownian::read_restart(FILE *fp) +{ + read_restart_settings(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(&cut_inner[i][j],sizeof(double),1,fp); + fread(&cut[i][j],sizeof(double),1,fp); + } + MPI_Bcast(&cut_inner[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairBrownian::write_restart_settings(FILE *fp) +{ + fwrite(&mu,sizeof(double),1,fp); + fwrite(&flaglog,sizeof(int),1,fp); + fwrite(&flagfld,sizeof(int),1,fp); + fwrite(&cut_inner_global,sizeof(double),1,fp); + fwrite(&cut_global,sizeof(double),1,fp); + fwrite(&t_target,sizeof(double),1,fp); + fwrite(&seed,sizeof(int),1,fp); + fwrite(&offset_flag,sizeof(int),1,fp); + fwrite(&mix_flag,sizeof(int),1,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairBrownian::read_restart_settings(FILE *fp) +{ + int me = comm->me; + if (me == 0) { + fread(&mu,sizeof(double),1,fp); + fread(&flaglog,sizeof(int),1,fp); + fread(&flagfld,sizeof(int),1,fp); + fread(&cut_inner_global,sizeof(double),1,fp); + fread(&cut_global,sizeof(double),1,fp); + fread(&t_target, sizeof(double),1,fp); + fread(&seed, sizeof(int),1,fp); + fread(&offset_flag,sizeof(int),1,fp); + fread(&mix_flag,sizeof(int),1,fp); + } + MPI_Bcast(&mu,1,MPI_DOUBLE,0,world); + MPI_Bcast(&flaglog,1,MPI_INT,0,world); + MPI_Bcast(&flagfld,1,MPI_INT,0,world); + MPI_Bcast(&cut_inner_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&t_target,1,MPI_DOUBLE,0,world); + MPI_Bcast(&seed,1,MPI_INT,0,world); + MPI_Bcast(&offset_flag,1,MPI_INT,0,world); + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); + + // additional setup based on restart parameters + + delete random; + random = new RanMars(lmp,seed + comm->me); +} + +/* ----------------------------------------------------------------------*/ + +void PairBrownian::set_3_orthogonal_vectors(double p1[3], + double p2[3], double p3[3]) +{ + double norm; + int ix,iy,iz; + + // find the index of maximum magnitude and store it in iz + + if (fabs(p1[0]) > fabs(p1[1])) { + iz=0; + ix=1; + iy=2; + } else { + iz=1; + ix=2; + iy=0; + } + + if (iz==0) { + if (fabs(p1[0]) < fabs(p1[2])) { + iz = 2; + ix = 0; + iy = 1; + } + } else { + if (fabs(p1[1]) < fabs(p1[2])) { + iz = 2; + ix = 0; + iy = 1; + } + } + + // set p2 arbitrarily such that it's orthogonal to p1 + + p2[ix]=1.0; + p2[iy]=1.0; + p2[iz] = -(p1[ix]*p2[ix] + p1[iy]*p2[iy])/p1[iz]; + + // normalize p2 + + norm = sqrt(pow(p2[0],2) + pow(p2[1],2) + pow(p2[2],2)); + + p2[0] = p2[0]/norm; + p2[1] = p2[1]/norm; + p2[2] = p2[2]/norm; + + // Set p3 by taking the cross product p3=p2xp1 + + p3[0] = p1[1]*p2[2] - p1[2]*p2[1]; + p3[1] = p1[2]*p2[0] - p1[0]*p2[2]; + p3[2] = p1[0]*p2[1] - p1[1]*p2[0]; +} diff --git a/src/FLD/pair_brownian.h b/src/FLD/pair_brownian.h new file mode 100644 index 0000000000..0ec34215fb --- /dev/null +++ b/src/FLD/pair_brownian.h @@ -0,0 +1,58 @@ +/* ---------------------------------------------------------------------- + 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(brownian,PairBrownian) + +#else + +#ifndef LMP_PAIR_BROWNIAN_H +#define LMP_PAIR_BROWNIAN_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairBrownian : public Pair { + public: + PairBrownian(class LAMMPS *); + virtual ~PairBrownian(); + virtual void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + virtual double init_one(int, int); + virtual void init_style(); + void write_restart(FILE *); + void read_restart(FILE *); + void write_restart_settings(FILE *); + void read_restart_settings(FILE *); + + protected: + double cut_inner_global,cut_global; + double t_target,mu; + int flaglog,flagfld; + int seed; + double **cut_inner,**cut; + double R0,RT0; + + class RanMars *random; + + void set_3_orthogonal_vectors(double*,double*,double*); + void allocate(); +}; + +} + +#endif +#endif diff --git a/src/FLD/pair_brownian_poly.cpp b/src/FLD/pair_brownian_poly.cpp new file mode 100644 index 0000000000..aec8775b67 --- /dev/null +++ b/src/FLD/pair_brownian_poly.cpp @@ -0,0 +1,338 @@ +/* ---------------------------------------------------------------------- + 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: Amit Kumar and Michael Bybee (UIUC) + Dave Heine (Corning), polydispersity +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_brownian_poly.h" +#include "atom.h" +#include "atom_vec.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "domain.h" +#include "update.h" +#include "memory.h" +#include "random_mars.h" +#include "math_const.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +/* ---------------------------------------------------------------------- */ + +PairBrownianPoly::PairBrownianPoly(LAMMPS *lmp) : PairBrownian(lmp) +{ + no_virial_fdotr_compute = 1; +} + +/* ---------------------------------------------------------------------- */ + +void PairBrownianPoly::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,tx,ty,tz; + double rsq,r,h_sep,beta0,beta1,radi,radj; + int *ilist,*jlist,*numneigh,**firstneigh; + + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + double **torque = atom->torque; + double *radius = atom->radius; + int *type = atom->type; + int nlocal = atom->nlocal; + + double vxmu2f = force->vxmu2f; + int overlaps = 0; + double randr; + double prethermostat; + double xl[3],a_sq,a_sh,a_pu,Fbmag; + double p1[3],p2[3],p3[3]; + + // scale factor for Brownian moments + + prethermostat = sqrt(24.0*force->boltz*t_target/update->dt); + prethermostat *= sqrt(force->vxmu2f/force->ftm2v/force->mvv2e); + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + radi = radius[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + // FLD contribution to force and torque due to isotropic terms + + if (flagfld) { + f[i][0] += prethermostat*sqrt(R0*radi)*(random->uniform()-0.5); + f[i][1] += prethermostat*sqrt(R0*radi)*(random->uniform()-0.5); + f[i][2] += prethermostat*sqrt(R0*radi)*(random->uniform()-0.5); + if (flaglog) { + torque[i][0] += prethermostat*sqrt(RT0*pow(radi,3)) * + (random->uniform()-0.5); + torque[i][1] += prethermostat*sqrt(RT0*pow(radi,3)) * + (random->uniform()-0.5); + torque[i][2] += prethermostat*sqrt(RT0*pow(radi,3)) * + (random->uniform()-0.5); + } + } + + 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; + jtype = type[j]; + radj = radius[j]; + + if (rsq < cutsq[itype][jtype]) { + r = sqrt(rsq); + + // scalar resistances a_sq and a_sh + + h_sep = r - radi-radj; + + // check for overlaps + + if (h_sep < 0.0) overlaps++; + + // if less than minimum gap, use minimum gap instead + + if (r < cut_inner[itype][jtype]) + h_sep = cut_inner[itype][jtype] - radi-radj; + + // scale h_sep by radi + + h_sep = h_sep/radi; + beta0 = radj/radi; + beta1 = 1.0 + beta0; + + // scalar resistances + + if (flaglog) { + a_sq = beta0*beta0/beta1/beta1/h_sep + + (1.0+7.0*beta0+beta0*beta0)/5.0/pow(beta1,3)*log(1.0/h_sep); + a_sq += (1.0+18.0*beta0-29.0*beta0*beta0+18.0*pow(beta0,3) + + pow(beta0,4))/21.0/pow(beta1,4)*h_sep*log(1.0/h_sep); + a_sq *= 6.0*MY_PI*mu*radi; + a_sh = 4.0*beta0*(2.0+beta0+2.0*beta0*beta0)/15.0/pow(beta1,3) * + log(1.0/h_sep); + a_sh += 4.0*(16.0-45.0*beta0+58.0*beta0*beta0-45.0*pow(beta0,3) + + 16.0*pow(beta0,4))/375.0/pow(beta1,4) * + h_sep*log(1.0/h_sep); + a_sh *= 6.0*MY_PI*mu*radi; + a_pu = beta0*(4.0+beta0)/10.0/beta1/beta1*log(1.0/h_sep); + a_pu += (32.0-33.0*beta0+83.0*beta0*beta0+43.0 * + pow(beta0,3))/250.0/pow(beta1,3)*h_sep*log(1.0/h_sep); + a_pu *= 8.0*MY_PI*mu*pow(radi,3); + + } else a_sq = 6.0*MY_PI*mu*radi*(beta0*beta0/beta1/beta1/h_sep); + + // generate the Pairwise Brownian Force: a_sq + + Fbmag = prethermostat*sqrt(a_sq); + + // generate a random number + + randr = random->uniform()-0.5; + + // contribution due to Brownian motion + + fx = Fbmag*randr*delx/r; + fy = Fbmag*randr*dely/r; + fz = Fbmag*randr*delz/r; + + // add terms due to a_sh + + if (flaglog) { + + // generate two orthogonal vectors to the line of centers + + p1[0] = delx/r; p1[1] = dely/r; p1[2] = delz/r; + set_3_orthogonal_vectors(p1,p2,p3); + + // magnitude + + Fbmag = prethermostat*sqrt(a_sh); + + // force in each of the two directions + + randr = random->uniform()-0.5; + fx += Fbmag*randr*p2[0]; + fy += Fbmag*randr*p2[1]; + fz += Fbmag*randr*p2[2]; + + randr = random->uniform()-0.5; + fx += Fbmag*randr*p3[0]; + fy += Fbmag*randr*p3[1]; + fz += Fbmag*randr*p3[2]; + } + + // scale forces to appropriate units + + fx = vxmu2f*fx; + fy = vxmu2f*fy; + fz = vxmu2f*fz; + + // sum to total Force + + f[i][0] -= fx; + f[i][1] -= fy; + f[i][2] -= fz; + + // torque due to the Brownian Force + + if (flaglog) { + + // location of the point of closest approach on I from its center + + xl[0] = -delx/r*radi; + xl[1] = -dely/r*radi; + xl[2] = -delz/r*radi; + + // torque = xl_cross_F + + tx = xl[1]*fz - xl[2]*fy; + ty = xl[2]*fx - xl[0]*fz; + tz = xl[0]*fy - xl[1]*fx; + + // torque is same on both particles + + torque[i][0] -= tx; + torque[i][1] -= ty; + torque[i][2] -= tz; + + // torque due to a_pu + + Fbmag = prethermostat*sqrt(a_pu); + + // force in each direction + + randr = random->uniform()-0.5; + tx = Fbmag*randr*p2[0]; + ty = Fbmag*randr*p2[1]; + tz = Fbmag*randr*p2[2]; + + randr = random->uniform()-0.5; + tx += Fbmag*randr*p3[0]; + ty += Fbmag*randr*p3[1]; + tz += Fbmag*randr*p3[2]; + + // torque has opposite sign on two particles + + torque[i][0] -= tx; + torque[i][1] -= ty; + torque[i][2] -= tz; + + } + + // set j = nlocal so that only I gets tallied + + if (evflag) ev_tally_xyz(i,nlocal,nlocal,0, + 0.0,0.0,-fx,-fy,-fz,delx,dely,delz); + } + } + } +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairBrownianPoly::init_style() +{ + if (force->newton_pair == 1) + error->all(FLERR,"Pair brownian/poly requires newton pair off"); + if (!atom->sphere_flag) + error->all(FLERR,"Pair brownian/poly requires atom style sphere"); + + // insure all particles are finite-size + // for pair hybrid, should limit test to types using the pair style + + double *radius = atom->radius; + int *type = atom->type; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (radius[i] == 0.0) + error->one(FLERR,"Pair brownian/poly requires extended particles"); + + int irequest = neighbor->request(this); + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->full = 1; + + // set the isotropic constants that depend on the volume fraction + // vol_T = total volume + + double vol_T = domain->xprd*domain->yprd*domain->zprd; + + // vol_P = volume of particles, assuming mono-dispersity + // vol_f = volume fraction + + double volP = 0.0; + + for (int i = 0; i < nlocal; i++) + volP += (4.0/3.0)*MY_PI*pow(atom->radius[i],3); + double vol_P; + MPI_Allreduce(&volP,&vol_P,1,MPI_DOUBLE,MPI_SUM,world); + double vol_f = vol_P/vol_T; + + // set isotropic constants + + if (flaglog == 0) { + R0 = 6*MY_PI*mu*(1.0 + 2.16*vol_f); + RT0 = 8*MY_PI*mu; + } else { + R0 = 6*MY_PI*mu*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f); + RT0 = 8*MY_PI*mu*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); + } +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairBrownianPoly::init_one(int i, int j) +{ + if (setflag[i][j] == 0) { + cut_inner[i][j] = mix_distance(cut_inner[i][i],cut_inner[j][j]); + cut[i][j] = mix_distance(cut[i][i],cut[j][j]); + } + + cut_inner[j][i] = cut_inner[i][j]; + return cut[i][j]; +} diff --git a/src/FLD/pair_brownian_poly.h b/src/FLD/pair_brownian_poly.h new file mode 100644 index 0000000000..e92d926df8 --- /dev/null +++ b/src/FLD/pair_brownian_poly.h @@ -0,0 +1,39 @@ +/* ---------------------------------------------------------------------- + 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(brownian/poly,PairBrownianPoly) + +#else + +#ifndef LMP_PAIR_BROWNIAN_POLY_H +#define LMP_PAIR_BROWNIAN_POLY_H + +#include "pair_brownian.h" + +namespace LAMMPS_NS { + +class PairBrownianPoly : public PairBrownian { + public: + PairBrownianPoly(class LAMMPS *); + ~PairBrownianPoly() {} + void compute(int, int); + double init_one(int, int); + void init_style(); +}; + +} + +#endif +#endif diff --git a/src/FLD/pair_lubricate.cpp b/src/FLD/pair_lubricate.cpp new file mode 100755 index 0000000000..904b8e1c17 --- /dev/null +++ b/src/FLD/pair_lubricate.cpp @@ -0,0 +1,723 @@ +/* ---------------------------------------------------------------------- + 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: Randy Schunk (SNL) + Amit Kumar and Michael Bybee (UIUC) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_lubricate.h" +#include "atom.h" +#include "atom_vec.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "domain.h" +#include "update.h" +#include "modify.h" +#include "fix.h" +#include "fix_deform.h" +#include "memory.h" +#include "random_mars.h" +#include "math_const.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +// same as fix_deform.cpp + +enum{NO_REMAP,X_REMAP,V_REMAP}; + +/* ---------------------------------------------------------------------- */ + +PairLubricate::PairLubricate(LAMMPS *lmp) : Pair(lmp) +{ + single_enable = 0; + + // set comm size needed by this Pair + + comm_forward = 6; +} + +/* ---------------------------------------------------------------------- */ + +PairLubricate::~PairLubricate() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + memory->destroy(cut); + memory->destroy(cut_inner); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairLubricate::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,fpair,fx,fy,fz,tx,ty,tz; + double rsq,r,h_sep,radi,tfmag; + double vr1,vr2,vr3,vnnr,vn1,vn2,vn3; + double vt1,vt2,vt3,wt1,wt2,wt3,wdotn; + double inertia,inv_inertia,vRS0; + double vi[3],vj[3],wi[3],wj[3],xl[3]; + double a_sq,a_sh,a_pu,Fbmag,del,delmin,eta; + int *ilist,*jlist,*numneigh,**firstneigh; + double lamda[3],vstream[3]; + + double vxmu2f = force->vxmu2f; + + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + double **omega = atom->omega; + double **angmom = atom->angmom; + double **torque = atom->torque; + double *radius = atom->radius; + double *mass = atom->mass; + double *rmass = atom->rmass; + int *type = atom->type; + int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; + + int overlaps = 0; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // subtract streaming component of velocity, omega, angmom + // assume fluid streaming velocity = box deformation rate + // vstream = (ux,uy,uz) + // ux = h_rate[0]*x + h_rate[5]*y + h_rate[4]*z + // uy = h_rate[1]*y + h_rate[3]*z + // uz = h_rate[2]*z + // omega_new = omega - curl(vstream)/2 + // angmom_new = angmom - I*curl(vstream)/2 + // Ef = (grad(vstream) + (grad(vstream))^T) / 2 + + if (shearing) { + double *h_rate = domain->h_rate; + double *h_ratelo = domain->h_ratelo; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + itype = type[i]; + radi = radius[i]; + + domain->x2lamda(x[i],lamda); + vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + + h_rate[4]*lamda[2] + h_ratelo[0]; + vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; + vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2]; + v[i][0] -= vstream[0]; + v[i][1] -= vstream[1]; + v[i][2] -= vstream[2]; + + omega[i][0] += 0.5*h_rate[3]; + omega[i][1] -= 0.5*h_rate[4]; + omega[i][2] += 0.5*h_rate[5]; + } + + // set Ef from h_rate in strain units + + Ef[0][0] = h_rate[0]/domain->xprd; + Ef[1][1] = h_rate[1]/domain->yprd; + Ef[2][2] = h_rate[2]/domain->zprd; + Ef[0][1] = Ef[1][0] = 0.5 * h_rate[5]/domain->yprd; + Ef[0][2] = Ef[2][0] = 0.5 * h_rate[4]/domain->zprd; + Ef[1][2] = Ef[2][1] = 0.5 * h_rate[3]/domain->zprd; + + // copy updated velocity/omega/angmom to the ghost particles + // no need to do this if not shearing since comm->ghost_velocity is set + + comm->forward_comm_pair(this); + } + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + radi = radius[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + // angular velocity + + wi[0] = omega[i][0]; + wi[1] = omega[i][1]; + wi[2] = omega[i][2]; + + // FLD contribution to force and torque due to isotropic terms + // FLD contribution to stress from isotropic RS0 + + if (flagfld) { + f[i][0] -= vxmu2f*R0*v[i][0]; + f[i][1] -= vxmu2f*R0*v[i][1]; + f[i][2] -= vxmu2f*R0*v[i][2]; + torque[i][0] -= vxmu2f*RT0*wi[0]; + torque[i][1] -= vxmu2f*RT0*wi[1]; + torque[i][2] -= vxmu2f*RT0*wi[2]; + + if (shearing && vflag_either) { + vRS0 = -vxmu2f * RS0; + v_tally_tensor(i,i,nlocal,newton_pair, + vRS0*Ef[0][0],vRS0*Ef[1][1],vRS0*Ef[2][2], + vRS0*Ef[0][1],vRS0*Ef[0][2],vRS0*Ef[1][2]); + } + } + + 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; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + r = sqrt(rsq); + + // angular momentum = I*omega = 2/5 * M*R^2 * omega + + wj[0] = omega[j][0]; + wj[1] = omega[j][1]; + wj[2] = omega[j][2]; + + // xl = point of closest approach on particle i from its center + + xl[0] = -delx/r*radi; + xl[1] = -dely/r*radi; + xl[2] = -delz/r*radi; + + // velocity at the point of closest approach on both particles + // v = v + omega_cross_xl - Ef.xl + + // particle i + + vi[0] = v[i][0] + (wi[1]*xl[2] - wi[2]*xl[1]) + - (Ef[0][0]*xl[0] + Ef[0][1]*xl[1] + Ef[0][2]*xl[2]); + + vi[1] = v[i][1] + (wi[2]*xl[0] - wi[0]*xl[2]) + - (Ef[1][0]*xl[0] + Ef[1][1]*xl[1] + Ef[1][2]*xl[2]); + + vi[2] = v[i][2] + (wi[0]*xl[1] - wi[1]*xl[0]) + - (Ef[2][0]*xl[0] + Ef[2][1]*xl[1] + Ef[2][2]*xl[2]); + + // particle j + + vj[0] = v[j][0] - (wj[1]*xl[2] - wj[2]*xl[1]) + + (Ef[0][0]*xl[0] + Ef[0][1]*xl[1] + Ef[0][2]*xl[2]); + + vj[1] = v[j][1] - (wj[2]*xl[0] - wj[0]*xl[2]) + + (Ef[1][0]*xl[0] + Ef[1][1]*xl[1] + Ef[1][2]*xl[2]); + + vj[2] = v[j][2] - (wj[0]*xl[1] - wj[1]*xl[0]) + + (Ef[2][0]*xl[0] + Ef[2][1]*xl[1] + Ef[2][2]*xl[2]); + + // scalar resistances XA and YA + + h_sep = r - 2.0*radi; + + // check for overlaps + + if (h_sep < 0.0) overlaps++; + + // if less than the minimum gap use the minimum gap instead + + if (r < cut_inner[itype][jtype]) + h_sep = cut_inner[itype][jtype] - 2.0*radi; + + // scale h_sep by radi + + h_sep = h_sep/radi; + + // scalar resistances + + if (flaglog) { + a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep + 9.0/40.0*log(1.0/h_sep)); + a_sh = 6.0*MY_PI*mu*radi*(1.0/6.0*log(1.0/h_sep)); + a_pu = 8.0*MY_PI*mu*pow(radi,3)*(3.0/160.0*log(1.0/h_sep)); + } else + a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep); + + // relative velocity at the point of closest approach + // includes fluid velocity + + vr1 = vi[0] - vj[0]; + vr2 = vi[1] - vj[1]; + vr3 = vi[2] - vj[2]; + + // normal component (vr.n)n + + vnnr = (vr1*delx + vr2*dely + vr3*delz)/r; + vn1 = vnnr*delx/r; + vn2 = vnnr*dely/r; + vn3 = vnnr*delz/r; + + // tangential component vr - (vr.n)n + + vt1 = vr1 - vn1; + vt2 = vr2 - vn2; + vt3 = vr3 - vn3; + + // force due to squeeze type motion + + fx = a_sq*vn1; + fy = a_sq*vn2; + fz = a_sq*vn3; + + // force due to all shear kind of motions + + if (flaglog) { + fx = fx + a_sh*vt1; + fy = fy + a_sh*vt2; + fz = fz + a_sh*vt3; + } + + // scale forces for appropriate units + + fx *= vxmu2f; + fy *= vxmu2f; + fz *= vxmu2f; + + // add to total force + + f[i][0] -= fx; + f[i][1] -= fy; + f[i][2] -= fz; + + if (newton_pair || j < nlocal) { + f[j][0] += fx; + f[j][1] += fy; + f[j][2] += fz; + } + + // torque due to this force + + if (flaglog) { + tx = xl[1]*fz - xl[2]*fy; + ty = xl[2]*fx - xl[0]*fz; + tz = xl[0]*fy - xl[1]*fx; + + torque[i][0] -= vxmu2f*tx; + torque[i][1] -= vxmu2f*ty; + torque[i][2] -= vxmu2f*tz; + + if (newton_pair || j < nlocal) { + torque[j][0] -= vxmu2f*tx; + torque[j][1] -= vxmu2f*ty; + torque[j][2] -= vxmu2f*tz; + } + + // torque due to a_pu + + wdotn = ((wi[0]-wj[0])*delx + (wi[1]-wj[1])*dely + + (wi[2]-wj[2])*delz)/r; + wt1 = (wi[0]-wj[0]) - wdotn*delx/r; + wt2 = (wi[1]-wj[1]) - wdotn*dely/r; + wt3 = (wi[2]-wj[2]) - wdotn*delz/r; + + tx = a_pu*wt1; + ty = a_pu*wt2; + tz = a_pu*wt3; + + torque[i][0] -= vxmu2f*tx; + torque[i][1] -= vxmu2f*ty; + torque[i][2] -= vxmu2f*tz; + + if (newton_pair || j < nlocal) { + torque[j][0] += vxmu2f*tx; + torque[j][1] += vxmu2f*ty; + torque[j][2] += vxmu2f*tz; + } + } + + if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair, + 0.0,0.0,-fx,-fy,-fz,delx,dely,delz); + } + } + } + + // restore streaming component of velocity, omega, angmom + + if (shearing) { + double *h_rate = domain->h_rate; + double *h_ratelo = domain->h_ratelo; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + itype = type[i]; + radi = radius[i]; + + domain->x2lamda(x[i],lamda); + vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + + h_rate[4]*lamda[2] + h_ratelo[0]; + vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; + vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2]; + v[i][0] += vstream[0]; + v[i][1] += vstream[1]; + v[i][2] += vstream[2]; + + omega[i][0] -= 0.5*h_rate[3]; + omega[i][1] += 0.5*h_rate[4]; + omega[i][2] -= 0.5*h_rate[5]; + } + } + + // to DEBUG: set print_overlaps to 1 + + int print_overlaps = 0; + if (print_overlaps) { + int overlaps_all; + MPI_Allreduce(&overlaps,&overlaps_all,1,MPI_INT,MPI_SUM,world); + if (overlaps_all && comm->me == 0) + printf("Number of overlaps = %d\n",overlaps); + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairLubricate::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(cut,n+1,n+1,"pair:cut"); + memory->create(cut_inner,n+1,n+1,"pair:cut_inner"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairLubricate::settings(int narg, char **arg) +{ + if (narg != 5) error->all(FLERR,"Illegal pair_style command"); + + mu = atof(arg[0]); + flaglog = atoi(arg[1]); + flagfld = atoi(arg[2]); + cut_inner_global = atof(arg[3]); + cut_global = atof(arg[4]); + + // reset cutoffs that have been explicitly set + + if (allocated) { + for (int i = 1; i <= atom->ntypes; i++) + for (int j = i+1; j <= atom->ntypes; j++) + if (setflag[i][j]) { + cut_inner[i][j] = cut_inner_global; + cut[i][j] = cut_global; + } + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairLubricate::coeff(int narg, char **arg) +{ + if (narg != 2 && narg != 4) + 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 cut_inner_one = cut_inner_global; + double cut_one = cut_global; + if (narg == 4) { + cut_inner_one = atof(arg[2]); + cut_one = atof(arg[3]); + } + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + cut_inner[i][j] = cut_inner_one; + cut[i][j] = cut_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairLubricate::init_style() +{ + if (!atom->sphere_flag) + error->all(FLERR,"Pair lubricate requires atom style sphere"); + if (comm->ghost_velocity == 0) + error->all(FLERR,"Pair lubricate requires ghost atoms store velocity"); + + int irequest = neighbor->request(this); + + // require that atom radii are identical within each type + // require monodisperse system with same radii for all types + + double rad,radtype; + for (int i = 1; i <= atom->ntypes; i++) { + if (!atom->radius_consistency(i,radtype)) + error->all(FLERR,"Pair lubricate requires monodisperse particles"); + if (i > 1 && radtype != rad) + error->all(FLERR,"Pair lubricate requires monodisperse particles"); + rad = radtype; + } + + // set the isotropic constants that depend on the volume fraction + // vol_T = total volume + + double vol_T = domain->xprd*domain->yprd*domain->zprd; + + // vol_P = volume of particles, assuming monodispersity + // vol_f = volume fraction + + double vol_P = atom->natoms*(4.0/3.0)*MY_PI*pow(rad,3); + double vol_f = vol_P/vol_T; + + // set isotropic constants for FLD + + if (flaglog == 0) { + R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f); + RT0 = 8*MY_PI*mu*pow(rad,3); + RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f); + } else { + R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f); + RT0 = 8*MY_PI*mu*pow(rad,3)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); + RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f); + } + + // check for fix deform, if exists it must use "remap v" + + shearing = 0; + for (int i = 0; i < modify->nfix; i++) + if (strcmp(modify->fix[i]->style,"deform") == 0) { + shearing = 1; + if (((FixDeform *) modify->fix[i])->remapflag != V_REMAP) + error->all(FLERR,"Using pair lubricate with inconsistent " + "fix deform remap option"); + } + + // set Ef = 0 since used whether shearing or not + + Ef[0][0] = Ef[0][1] = Ef[0][2] = 0.0; + Ef[1][0] = Ef[1][1] = Ef[1][2] = 0.0; + Ef[2][0] = Ef[2][1] = Ef[2][2] = 0.0; +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairLubricate::init_one(int i, int j) +{ + if (setflag[i][j] == 0) { + cut_inner[i][j] = mix_distance(cut_inner[i][i],cut_inner[j][j]); + cut[i][j] = mix_distance(cut[i][i],cut[j][j]); + } + + cut_inner[j][i] = cut_inner[i][j]; + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLubricate::write_restart(FILE *fp) +{ + write_restart_settings(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(&cut_inner[i][j],sizeof(double),1,fp); + fwrite(&cut[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLubricate::read_restart(FILE *fp) +{ + read_restart_settings(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(&cut_inner[i][j],sizeof(double),1,fp); + fread(&cut[i][j],sizeof(double),1,fp); + } + MPI_Bcast(&cut_inner[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLubricate::write_restart_settings(FILE *fp) +{ + fwrite(&mu,sizeof(double),1,fp); + fwrite(&flaglog,sizeof(int),1,fp); + fwrite(&flagfld,sizeof(int),1,fp); + fwrite(&cut_inner_global,sizeof(double),1,fp); + fwrite(&cut_global,sizeof(double),1,fp); + fwrite(&offset_flag,sizeof(int),1,fp); + fwrite(&mix_flag,sizeof(int),1,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLubricate::read_restart_settings(FILE *fp) +{ + int me = comm->me; + if (me == 0) { + fread(&mu,sizeof(double),1,fp); + fread(&flaglog,sizeof(int),1,fp); + fread(&flagfld,sizeof(int),1,fp); + fread(&cut_inner_global,sizeof(double),1,fp); + fread(&cut_global,sizeof(double),1,fp); + fread(&offset_flag,sizeof(int),1,fp); + fread(&mix_flag,sizeof(int),1,fp); + } + MPI_Bcast(&mu,1,MPI_DOUBLE,0,world); + MPI_Bcast(&flaglog,1,MPI_INT,0,world); + MPI_Bcast(&flagfld,1,MPI_INT,0,world); + MPI_Bcast(&cut_inner_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&offset_flag,1,MPI_INT,0,world); + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); +} + +/* ---------------------------------------------------------------------- */ + +int PairLubricate::pack_comm(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + + double **v = atom->v; + double **omega = atom->omega; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + buf[m++] = omega[j][0]; + buf[m++] = omega[j][1]; + buf[m++] = omega[j][2]; + } + + return 6; +} + +/* ---------------------------------------------------------------------- */ + +void PairLubricate::unpack_comm(int n, int first, double *buf) +{ + int i,m,last; + + double **v = atom->v; + double **omega = atom->omega; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + omega[i][0] = buf[m++]; + omega[i][1] = buf[m++]; + omega[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- + check if name is recognized, return integer index for that name + if name not recognized, return -1 + if type pair setting, return -2 if no type pairs are set +------------------------------------------------------------------------- */ + +int PairLubricate::pre_adapt(char *name, int ilo, int ihi, int jlo, int jhi) +{ + if (strcmp(name,"mu") == 0) return 0; + return -1; +} + +/* ---------------------------------------------------------------------- + adapt parameter indexed by which + change all pair variables affected by the reset parameter + if type pair setting, set I-J and J-I coeffs +------------------------------------------------------------------------- */ + +void PairLubricate::adapt(int which, int ilo, int ihi, int jlo, int jhi, + double value) +{ + mu = value; +} diff --git a/src/FLD/pair_lubricate.h b/src/FLD/pair_lubricate.h new file mode 100644 index 0000000000..10cd0c7ecf --- /dev/null +++ b/src/FLD/pair_lubricate.h @@ -0,0 +1,59 @@ +/* ---------------------------------------------------------------------- + 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(lubricate,PairLubricate) + +#else + +#ifndef LMP_PAIR_LUBRICATE_H +#define LMP_PAIR_LUBRICATE_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairLubricate : public Pair { + public: + PairLubricate(class LAMMPS *); + virtual ~PairLubricate(); + virtual void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + double init_one(int, int); + virtual void init_style(); + void write_restart(FILE *); + void read_restart(FILE *); + void write_restart_settings(FILE *); + void read_restart_settings(FILE *); + int pre_adapt(char *, int, int, int, int); + void adapt(int, int, int, int, int, double); + + int pack_comm(int, int *, double *, int, int *); + void unpack_comm(int, int, double *); + + protected: + double mu,cut_inner_global,cut_global; + int flaglog,flagfld,shearing; + double Ef[3][3]; + double R0,RT0,RS0; + double **cut_inner,**cut; + + void allocate(); +}; + +} + +#endif +#endif diff --git a/src/FLD/pair_lubricateU.cpp b/src/FLD/pair_lubricateU.cpp new file mode 100644 index 0000000000..9cde9568a7 --- /dev/null +++ b/src/FLD/pair_lubricateU.cpp @@ -0,0 +1,1933 @@ +/* ---------------------------------------------------------------------- + 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: Amit Kumar and Michael Bybee (UIUC) +------------------------------------------------------------------------- */ + +#include "mpi.h" +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_lubricateU.h" +#include "atom.h" +#include "atom_vec.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "domain.h" +#include "update.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +#define TOL 1E-4 // tolerance for conjugate gradient + +/* ---------------------------------------------------------------------- */ + +PairLubricateU::PairLubricateU(LAMMPS *lmp) : Pair(lmp) +{ + single_enable = 0; + + // pair lubricateU cannot compute virial as F dot r + // due to how drag forces are applied to atoms + // correct method is how per-atom virial does it + + no_virial_fdotr_compute = 1; + + nmax = 0; + fl = Tl = xl = NULL; + + cgmax = 0; + bcg = xcg = rcg = rcg1 = pcg = RU = NULL; + + // set comm size needed by this Pair + + comm_forward = 6; +} + +/* ---------------------------------------------------------------------- */ + +PairLubricateU::~PairLubricateU() +{ + memory->destroy(fl); + memory->destroy(Tl); + memory->destroy(xl); + + memory->destroy(bcg); + memory->destroy(xcg); + memory->destroy(rcg); + memory->destroy(rcg1); + memory->destroy(pcg); + memory->destroy(RU); + + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + memory->destroy(cut); + memory->destroy(cut_inner); + } +} + +/* ---------------------------------------------------------------------- + It first has to solve for the velocity of the particles such that + the net force on the particles is zero. NOTE: it has to be the last + type of pair interaction specified in the input file. Also, it + assumes that no other types of interactions, like k-space, is + present. As already mentioned, the net force on the particles after + this pair interaction would be identically zero. + ---------------------------------------------------------------------- */ + +void PairLubricateU::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + double **omega = atom->omega; + double **torque = atom->torque; + int *type = atom->type; + int nlocal = atom->nlocal; + int nghost = atom->nghost; + int nall = nlocal + nghost; + + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + // skip compute() if called from integrate::setup() + // this is b/c do not want compute() to update velocities twice on a restart + // when restarting, call compute on step N (last step of prev run), + // again on step N (setup of restart run), + // then on step N+1 (first step of restart) + // so this is one extra time which leads to bad dynamics + + if (update->setupflag) return; + + // grow per-atom arrays if necessary + // need to be atom->nmax in length + + if (atom->nmax > nmax) { + memory->destroy(fl); + memory->destroy(Tl); + memory->destroy(xl); + nmax = atom->nmax; + memory->create(fl,nmax,3,"pair:fl"); + memory->create(Tl,nmax,3,"pair:Tl"); + memory->create(xl,nmax,3,"pair:xl"); + } + + // Added to implement midpoint integration scheme + // Save force, torque found so far. Also save the positions + + for (i=0;ix; + double **v = atom->v; + double **f = atom->f; + double **omega = atom->omega; + double **torque = atom->torque; + double *radius = atom->radius; + double *rmass = atom->rmass; + int *type = atom->type; + int nlocal = atom->nlocal; + + int newton_pair = force->newton_pair; + double vxmu2f = force->vxmu2f; + double inv_inertia,mo_inertia; + int *ilist; + + double radi; + int nprocs = comm->nprocs; + + inum = list->inum; + ilist = list->ilist; + + if (6*inum > cgmax) { + memory->destroy(bcg); + memory->destroy(xcg); + memory->destroy(rcg); + memory->destroy(rcg1); + memory->destroy(pcg); + memory->destroy(RU); + cgmax = 6*inum; + memory->create(bcg,cgmax,"pair:bcg"); + memory->create(xcg,cgmax,"pair:bcg"); + memory->create(rcg,cgmax,"pair:bcg"); + memory->create(rcg1,cgmax,"pair:bcg"); + memory->create(pcg,cgmax,"pair:bcg"); + memory->create(RU,cgmax,"pair:bcg"); + } + + double alpha,beta; + double normi,error,normig; + double send[2],recv[2],rcg_dot_rcg; + + // First compute R_FE*E + + compute_RE(); + + // Reverse communication of forces and torques to + // accumulate the net force on each of the particles + + if (newton_pair) comm->reverse_comm(); + + // CONJUGATE GRADIENT + // Find the right hand side= -ve of all forces/torques + // b = 6*Npart in overall size + + for(ii = 0; ii < inum; ii++) { + i = ilist[ii]; + for (j = 0; j < 3; j++) { + bcg[6*ii+j] = -f[i][j]; + bcg[6*ii+j+3] = -torque[i][j]; + } + } + + // Start solving the equation : F^H = -F^P -F^B - F^H_{Ef} + // Store initial guess for velocity and angular-velocities/angular momentum + // NOTE velocities and angular velocities are assumed relative to the fluid + + for (ii=0;iiforward_comm_pair(this); + + // Find initial residual + + compute_RU(); + + // reverse communication of forces and torques + + if (newton_pair) comm->reverse_comm(); + + copy_uo_vec(inum,f,torque,RU); + + for (i=0;i<6*inum;i++) + rcg[i] = bcg[i] - RU[i]; + + // Set initial conjugate direction + + for (i=0;i<6*inum;i++) + pcg[i] = rcg[i]; + + // Find initial norm of the residual or norm of the RHS (either is fine) + + normi = dot_vec_vec(6*inum,bcg,bcg); + + MPI_Allreduce(&normi,&normig,1,MPI_DOUBLE,MPI_SUM,world); + + // Loop until convergence + + do { + // find R*p + + copy_vec_uo(inum,pcg,v,omega); + + // set velocities for ghost particles + + comm->forward_comm_pair(this); + + compute_RU(); + + // reverse communication of forces and torques + + if (newton_pair) comm->reverse_comm(); + + + copy_uo_vec(inum,f,torque,RU); + + // Find alpha + + send[0] = dot_vec_vec(6*inum,rcg,rcg); + send[1] = dot_vec_vec(6*inum,RU,pcg); + + MPI_Allreduce(send,recv,2,MPI_DOUBLE,MPI_SUM,world); + + alpha = recv[0]/recv[1]; + rcg_dot_rcg = recv[0]; + + // Find new x + + for (i=0;i<6*inum;i++) + xcg[i] = xcg[i] + alpha*pcg[i]; + + // find new residual + + for (i=0;i<6*inum;i++) + rcg1[i] = rcg[i] - alpha*RU[i]; + + // find beta + + send[0] = dot_vec_vec(6*inum,rcg1,rcg1); + + MPI_Allreduce(send,recv,1,MPI_DOUBLE,MPI_SUM,world); + + beta = recv[0]/rcg_dot_rcg; + + // Find new conjugate direction + + for (i=0;i<6*inum;i++) + pcg[i] = rcg1[i] + beta*pcg[i]; + + for (i=0;i<6*inum;i++) + rcg[i] = rcg1[i]; + + // Find relative error + + error = sqrt(recv[0]/normig); + + } while (error > TOL); + + // update the final converged velocities in respective arrays + + copy_vec_uo(inum,xcg,v,omega); + + // set velocities for ghost particles + + comm->forward_comm_pair(this); + + // Find actual particle's velocities from relative velocities + // Only non-zero component of fluid's vel : vx=gdot*y and wz=-gdot/2 + + for (ii=0;iix; + double **v = atom->v; + double dtv = update->dt; + + for (i=0;iv; + double **f = atom->f; + double **omega = atom->omega; + double **torque = atom->torque; + double *radius = atom->radius; + double *rmass = atom->rmass; + int *type = atom->type; + int nlocal = atom->nlocal; + + int newton_pair = force->newton_pair; + double vxmu2f = force->vxmu2f; + double inv_inertia,mo_inertia; + int *ilist; + + double radi; + int nprocs = comm->nprocs; + + inum = list->inum; + ilist = list->ilist; + + double alpha,beta; + double normi,error,normig; + double send[2],recv[2],rcg_dot_rcg; + + // First compute R_FE*E + + compute_RE(x); + + // Reverse communication of forces and torques to + // accumulate the net force on each of the particles + + if (newton_pair) comm->reverse_comm(); + + // CONJUGATE GRADIENT + // Find the right hand side= -ve of all forces/torques + // b = 6*Npart in overall size + + for(ii = 0; ii < inum; ii++) { + i = ilist[ii]; + for (j = 0; j < 3; j++) { + bcg[6*ii+j] = -f[i][j]; + bcg[6*ii+j+3] = -torque[i][j]; + } + } + + // Start solving the equation : F^H = -F^P -F^B - F^H_{Ef} + // Store initial guess for velocity and angular-velocities/angular momentum + // NOTE velocities and angular velocities are assumed relative to the fluid + + for (ii=0;iiforward_comm_pair(this); + + // Find initial residual + + compute_RU(x); + + // reverse communication of forces and torques + + if (newton_pair) comm->reverse_comm(); + + copy_uo_vec(inum,f,torque,RU); + + for (i=0;i<6*inum;i++) + rcg[i] = bcg[i] - RU[i]; + + // Set initial conjugate direction + + for (i=0;i<6*inum;i++) + pcg[i] = rcg[i]; + + // Find initial norm of the residual or norm of the RHS (either is fine) + + normi = dot_vec_vec(6*inum,bcg,bcg); + + MPI_Allreduce(&normi,&normig,1,MPI_DOUBLE,MPI_SUM,world); + + // Loop until convergence + + do { + // find R*p + + copy_vec_uo(inum,pcg,v,omega); + + // set velocities for ghost particles + + comm->forward_comm_pair(this); + + compute_RU(x); + + // reverse communication of forces and torques + + if (newton_pair) comm->reverse_comm(); + + copy_uo_vec(inum,f,torque,RU); + + // Find alpha + + send[0] = dot_vec_vec(6*inum,rcg,rcg); + send[1] = dot_vec_vec(6*inum,RU,pcg); + + MPI_Allreduce(send,recv,2,MPI_DOUBLE,MPI_SUM,world); + + alpha = recv[0]/recv[1]; + rcg_dot_rcg = recv[0]; + + // Find new x + + for (i=0;i<6*inum;i++) + xcg[i] = xcg[i] + alpha*pcg[i]; + + // find new residual + + for (i=0;i<6*inum;i++) + rcg1[i] = rcg[i] - alpha*RU[i]; + + // find beta + + send[0] = dot_vec_vec(6*inum,rcg1,rcg1); + + MPI_Allreduce(send,recv,1,MPI_DOUBLE,MPI_SUM,world); + + beta = recv[0]/rcg_dot_rcg; + + // Find new conjugate direction + + for (i=0;i<6*inum;i++) + pcg[i] = rcg1[i] + beta*pcg[i]; + + for (i=0;i<6*inum;i++) + rcg[i] = rcg1[i]; + + // Find relative error + + error = sqrt(recv[0]/normig); + + } while (error > TOL); + + + // update the final converged velocities in respective arrays + + copy_vec_uo(inum,xcg,v,omega); + + // set velocities for ghost particles + + comm->forward_comm_pair(this); + + // Compute the viscosity/pressure + + if (evflag) compute_Fh(x); + + // Find actual particle's velocities from relative velocities + // Only non-zero component of fluid's vel : vx=gdot*y and wz=-gdot/2 + + for (ii=0;iiv; + double **f = atom->f; + double **omega = atom->omega; + double **torque = atom->torque; + double *radius = atom->radius; + int *type = atom->type; + int nlocal = atom->nlocal; + int nghost = atom->nghost; + int newton_pair = force->newton_pair; + + double vxmu2f = force->vxmu2f; + int overlaps = 0; + double vi[3],vj[3],wi[3],wj[3],xl[3],a_sq,a_sh,a_pu,Fbmag,del,delmin,eta; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // Set force to zero which is the final value after this pair interaction + + for (i=0;ireverse_comm(); // not really needed + + // Find additional contribution from the stresslets + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + radi = radius[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + // Find the contribution to stress from isotropic RS0 + // Set psuedo force to obtain the required contribution + // need to set delx and fy only + + fx = 0.0; delx = radi; + fy = vxmu2f*RS0*gdot/2.0/radi; dely = 0.0; + fz = 0.0; delz = 0.0; + if (evflag) + ev_tally_xyz(i,i,nlocal,newton_pair,0.0,0.0,-fx,-fy,-fz,delx,dely,delz); + + // Find angular velocity + + wi[0] = omega[i][0]; + wi[1] = omega[i][1]; + wi[2] = omega[i][2]; + + 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; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + r = sqrt(rsq); + + // Use omega directly if it exists, else angmom + // angular momentum = I*omega = 2/5 * M*R^2 * omega + + wj[0] = omega[j][0]; + wj[1] = omega[j][1]; + wj[2] = omega[j][2]; + + // loc of the point of closest approach on particle i from its cente + + xl[0] = -delx/r*radi; + xl[1] = -dely/r*radi; + xl[2] = -delz/r*radi; + + // velocity at the point of closest approach on both particles + // v = v + omega_cross_xl + + // particle i + + vi[0] = v[i][0] + (wi[1]*xl[2] - wi[2]*xl[1]); + vi[1] = v[i][1] + (wi[2]*xl[0] - wi[0]*xl[2]); + vi[2] = v[i][2] + (wi[0]*xl[1] - wi[1]*xl[0]); + + // particle j + + vj[0] = v[j][0] - (wj[1]*xl[2] - wj[2]*xl[1]); + vj[1] = v[j][1] - (wj[2]*xl[0] - wj[0]*xl[2]); + vj[2] = v[j][2] - (wj[0]*xl[1] - wj[1]*xl[0]); + + + // Relative velocity at the point of closest approach + // include contribution from Einf of the fluid + + vr1 = vi[0] - vj[0] - + 2.0*(Ef[0][0]*xl[0] + Ef[0][1]*xl[1] + Ef[0][2]*xl[2]); + vr2 = vi[1] - vj[1] - + 2.0*(Ef[1][0]*xl[0] + Ef[1][1]*xl[1] + Ef[1][2]*xl[2]); + vr3 = vi[2] - vj[2] - + 2.0*(Ef[2][0]*xl[0] + Ef[2][1]*xl[1] + Ef[2][2]*xl[2]); + + // Normal component (vr.n)n + + vnnr = (vr1*delx + vr2*dely + vr3*delz)/r; + vn1 = vnnr*delx/r; + vn2 = vnnr*dely/r; + vn3 = vnnr*delz/r; + + // Tangential component vr - (vr.n)n + + vt1 = vr1 - vn1; + vt2 = vr2 - vn2; + vt3 = vr3 - vn3; + + // Find the scalar resistances a_sq, a_sh and a_pu + + h_sep = r - 2.0*radi; + + // check for overlaps + + if (h_sep < 0.0) overlaps++; + + // If less than the minimum gap use the minimum gap instead + + if (r < cut_inner[itype][jtype]) + h_sep = cut_inner[itype][jtype] - 2.0*radi; + + // Scale h_sep by radi + + h_sep = h_sep/radi; + + // Scalar resistances + + if (flaglog) { + a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep + 9.0/40.0*log(1.0/h_sep)); + a_sh = 6.0*MY_PI*mu*radi*(1.0/6.0*log(1.0/h_sep)); + } else + a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep); + + // Find force due to squeeze type motion + + fx = a_sq*vn1; + fy = a_sq*vn2; + fz = a_sq*vn3; + + // Find force due to all shear kind of motions + + if (flaglog) { + fx = fx + a_sh*vt1; + fy = fy + a_sh*vt2; + fz = fz + a_sh*vt3; + } + + // Scale forces to obtain in appropriate units + + fx = vxmu2f*fx; + fy = vxmu2f*fy; + fz = vxmu2f*fz; + + if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair, + 0.0,0.0,-fx,-fy,-fz,delx,dely,delz); + } + } + } +} + +/* ---------------------------------------------------------------------- + computes R_FU * U +---------------------------------------------------------------------- */ + +void PairLubricateU::compute_RU() +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,fpair,fx,fy,fz,tx,ty,tz; + double rsq,r,h_sep,radi,tfmag; + double vr1,vr2,vr3,vnnr,vn1,vn2,vn3; + double vt1,vt2,vt3,wdotn,wt1,wt2,wt3; + double inv_inertia; + int *ilist,*jlist,*numneigh,**firstneigh; + + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + double **omega = atom->omega; + double **torque = atom->torque; + double *radius = atom->radius; + double *rmass = atom->rmass; + int *type = atom->type; + int nlocal = atom->nlocal; + int nghost = atom->nghost; + int newton_pair = force->newton_pair; + + double vxmu2f = force->vxmu2f; + int overlaps = 0; + double vi[3],vj[3],wi[3],wj[3],xl[3],a_sq,a_sh,a_pu,Fbmag,del,delmin,eta; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // Initialize f to zero + + for (i=0;iv; + double **f = atom->f; + double **omega = atom->omega; + double **torque = atom->torque; + double *radius = atom->radius; + int *type = atom->type; + int nlocal = atom->nlocal; + int nghost = atom->nghost; + int newton_pair = force->newton_pair; + + double vxmu2f = force->vxmu2f; + int overlaps = 0; + double vi[3],vj[3],wi[3],wj[3],xl[3],a_sq,a_sh,a_pu,Fbmag,del,delmin,eta; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // Initialize f to zero + + for (i=0;ix; + double **v = atom->v; + double **f = atom->f; + double **omega = atom->omega; + double **torque = atom->torque; + double *radius = atom->radius; + int *type = atom->type; + int nlocal = atom->nlocal; + int nghost = atom->nghost; + int newton_pair = force->newton_pair; + + double vxmu2f = force->vxmu2f; + int overlaps = 0; + double vi[3],vj[3],wi[3],wj[3],xl[3],a_sq,a_sh,a_pu,Fbmag,del,delmin,eta; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + radi = radius[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + // No contribution from isotropic terms due to E + + 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; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + r = sqrt(rsq); + + // loc of the point of closest approach on particle i from its center + + xl[0] = -delx/r*radi; + xl[1] = -dely/r*radi; + xl[2] = -delz/r*radi; + + // Find the scalar resistances a_sq and a_sh + + h_sep = r - 2.0*radi; + + // check for overlaps + + if(h_sep < 0.0) overlaps++; + + // If less than the minimum gap use the minimum gap instead + + if (r < cut_inner[itype][jtype]) + h_sep = cut_inner[itype][jtype] - 2.0*radi; + + // Scale h_sep by radi + + h_sep = h_sep/radi; + + // Scalar resistance for Squeeze type motions + + if (flaglog) + a_sq = 6*MY_PI*mu*radi*(1.0/4.0/h_sep + 9.0/40.0*log(1/h_sep)); + else + a_sq = 6*MY_PI*mu*radi*(1.0/4.0/h_sep); + + // Scalar resistance for Shear type motions + + if (flaglog) { + a_sh = 6*MY_PI*mu*radi*(1.0/6.0*log(1/h_sep)); + a_pu = 8.0*MY_PI*mu*pow(radi,3)*(3.0/160.0*log(1.0/h_sep)); + } + + // Relative velocity at the point of closest approach due to Ef only + + vr1 = -2.0*(Ef[0][0]*xl[0] + Ef[0][1]*xl[1] + Ef[0][2]*xl[2]); + vr2 = -2.0*(Ef[1][0]*xl[0] + Ef[1][1]*xl[1] + Ef[1][2]*xl[2]); + vr3 = -2.0*(Ef[2][0]*xl[0] + Ef[2][1]*xl[1] + Ef[2][2]*xl[2]); + + // Normal component (vr.n)n + + vnnr = (vr1*delx + vr2*dely + vr3*delz)/r; + vn1 = vnnr*delx/r; + vn2 = vnnr*dely/r; + vn3 = vnnr*delz/r; + + // Tangential component vr - (vr.n)n + + vt1 = vr1 - vn1; + vt2 = vr2 - vn2; + vt3 = vr3 - vn3; + + // Find force due to squeeze type motion + + fx = a_sq*vn1; + fy = a_sq*vn2; + fz = a_sq*vn3; + + // Find force due to all shear kind of motions + + if (flaglog) { + fx = fx + a_sh*vt1; + fy = fy + a_sh*vt2; + fz = fz + a_sh*vt3; + } + + // Scale forces to obtain in appropriate units + + fx = vxmu2f*fx; + fy = vxmu2f*fy; + fz = vxmu2f*fz; + + // Add to the total forc + + f[i][0] -= fx; + f[i][1] -= fy; + f[i][2] -= fz; + + if (newton_pair || j < nlocal) { + f[j][0] += fx; + f[j][1] += fy; + f[j][2] += fz; + } + + // Find torque due to this force + + if (flaglog) { + tx = xl[1]*fz - xl[2]*fy; + ty = xl[2]*fx - xl[0]*fz; + tz = xl[0]*fy - xl[1]*fx; + + // Why a scale factor ? + + torque[i][0] -= vxmu2f*tx; + torque[i][1] -= vxmu2f*ty; + torque[i][2] -= vxmu2f*tz; + + if (newton_pair || j < nlocal) { + torque[j][0] -= vxmu2f*tx; + torque[j][1] -= vxmu2f*ty; + torque[j][2] -= vxmu2f*tz; + } + + // NOTE No a_pu term needed as they add up to zero + } + } + } + } + + int print_overlaps = 0; + if (print_overlaps && overlaps) + printf("Number of overlaps=%d\n",overlaps); +} + +/* ---------------------------------------------------------------------- + This computes R_{FE}*E , where E is the rate of strain of tensor which is + known apriori, as it depends only on the known fluid velocity. + So, this part of the hydrodynamic interaction can be pre computed and + transferred to the RHS + ---------------------------------------------------------------------- */ + +void PairLubricateU::compute_RE(double **x) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,fpair,fx,fy,fz,tx,ty,tz; + double rsq,r,h_sep,radi,tfmag; + double vr1,vr2,vr3,vnnr,vn1,vn2,vn3; + double vt1,vt2,vt3; + double inv_inertia; + int *ilist,*jlist,*numneigh,**firstneigh; + + double **v = atom->v; + double **f = atom->f; + double **omega = atom->omega; + double **torque = atom->torque; + double *radius = atom->radius; + int *type = atom->type; + int nlocal = atom->nlocal; + int nghost = atom->nghost; + int newton_pair = force->newton_pair; + + double vxmu2f = force->vxmu2f; + int overlaps = 0; + double vi[3],vj[3],wi[3],wj[3],xl[3],a_sq,a_sh,a_pu,Fbmag,del,delmin,eta; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + radi = radius[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + // No contribution from isotropic terms due to E + + 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; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + r = sqrt(rsq); + + // loc of the point of closest approach on particle i from its center + + xl[0] = -delx/r*radi; + xl[1] = -dely/r*radi; + xl[2] = -delz/r*radi; + + // Find the scalar resistances a_sq and a_sh + + h_sep = r - 2.0*radi; + + // check for overlaps + + if(h_sep < 0.0) overlaps++; + + // If less than the minimum gap use the minimum gap instead + + if (r < cut_inner[itype][jtype]) + h_sep = cut_inner[itype][jtype] - 2.0*radi; + + // Scale h_sep by radi + + h_sep = h_sep/radi; + + // Scalar resistance for Squeeze type motions + + if (flaglog) + a_sq = 6*MY_PI*mu*radi*(1.0/4.0/h_sep + 9.0/40.0*log(1/h_sep)); + else + a_sq = 6*MY_PI*mu*radi*(1.0/4.0/h_sep); + + // Scalar resistance for Shear type motions + + if (flaglog) { + a_sh = 6*MY_PI*mu*radi*(1.0/6.0*log(1/h_sep)); + a_pu = 8.0*MY_PI*mu*pow(radi,3)*(3.0/160.0*log(1.0/h_sep)); + } + + // Relative velocity at the point of closest approach due to Ef only + + vr1 = -2.0*(Ef[0][0]*xl[0] + Ef[0][1]*xl[1] + Ef[0][2]*xl[2]); + vr2 = -2.0*(Ef[1][0]*xl[0] + Ef[1][1]*xl[1] + Ef[1][2]*xl[2]); + vr3 = -2.0*(Ef[2][0]*xl[0] + Ef[2][1]*xl[1] + Ef[2][2]*xl[2]); + + // Normal component (vr.n)n + + vnnr = (vr1*delx + vr2*dely + vr3*delz)/r; + vn1 = vnnr*delx/r; + vn2 = vnnr*dely/r; + vn3 = vnnr*delz/r; + + // Tangential component vr - (vr.n)n + + vt1 = vr1 - vn1; + vt2 = vr2 - vn2; + vt3 = vr3 - vn3; + + // Find force due to squeeze type motion + + fx = a_sq*vn1; + fy = a_sq*vn2; + fz = a_sq*vn3; + + // Find force due to all shear kind of motions + + if (flaglog) { + fx = fx + a_sh*vt1; + fy = fy + a_sh*vt2; + fz = fz + a_sh*vt3; + } + + // Scale forces to obtain in appropriate units + + fx = vxmu2f*fx; + fy = vxmu2f*fy; + fz = vxmu2f*fz; + + // Add to the total forc + + f[i][0] -= fx; + f[i][1] -= fy; + f[i][2] -= fz; + + if (newton_pair || j < nlocal) { + f[j][0] += fx; + f[j][1] += fy; + f[j][2] += fz; + } + + // Find torque due to this force + + if (flaglog) { + tx = xl[1]*fz - xl[2]*fy; + ty = xl[2]*fx - xl[0]*fz; + tz = xl[0]*fy - xl[1]*fx; + + // Why a scale factor ? + + torque[i][0] -= vxmu2f*tx; + torque[i][1] -= vxmu2f*ty; + torque[i][2] -= vxmu2f*tz; + + if (newton_pair || j < nlocal) { + torque[j][0] -= vxmu2f*tx; + torque[j][1] -= vxmu2f*ty; + torque[j][2] -= vxmu2f*tz; + } + + // NOTE No a_pu term needed as they add up to zero + } + } + } + } + + int print_overlaps = 0; + if (print_overlaps && overlaps) + printf("Number of overlaps=%d\n",overlaps); +} + + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairLubricateU::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + setflag = 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; + + cutsq = memory->create(cutsq,n+1,n+1,"pair:cutsq"); + + memory->create(cut,n+1,n+1,"pair:cut"); + memory->create(cut_inner,n+1,n+1,"pair:cut_inner"); +} + +/*----------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairLubricateU::settings(int narg, char **arg) +{ + if (narg != 5) error->all(FLERR,"Illegal pair_style command"); + + mu = atof(arg[0]); + flaglog = atoi(arg[1]); + cut_inner_global = atof(arg[2]); + cut_global = atof(arg[3]); + gdot = atof(arg[4]); + + // reset cutoffs that have been explicitly set + + if (allocated) { + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i+1; j <= atom->ntypes; j++) + if (setflag[i][j]) { + cut_inner[i][j] = cut_inner_global; + cut[i][j] = cut_global; + } + } + + // store the rate of strain tensor + + Ef[0][0] = 0.0; + Ef[0][1] = 0.5*gdot; + Ef[0][2] = 0.0; + Ef[1][0] = 0.5*gdot; + Ef[1][1] = 0.0; + Ef[1][2] = 0.0; + Ef[2][0] = 0.0; + Ef[2][1] = 0.0; + Ef[2][2] = 0.0; +} + +/*----------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairLubricateU::coeff(int narg, char **arg) +{ + if (narg != 2 && narg != 4) + 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 cut_inner_one = cut_inner_global; + double cut_one = cut_global; + if (narg == 4) { + cut_inner_one = atof(arg[2]); + cut_one = atof(arg[3]); + } + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + cut_inner[i][j] = cut_inner_one; + cut[i][j] = cut_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairLubricateU::init_style() +{ + if (!atom->sphere_flag) + error->all(FLERR,"Pair lubricateU requires atom style sphere"); + if (comm->ghost_velocity == 0) + error->all(FLERR,"Pair lubricateU requires ghost atoms store velocity"); + + int irequest = neighbor->request(this); + + // require that atom radii are identical within each type + // require monodisperse system with same radii for all types + + double rad,radtype; + for (int i = 1; i <= atom->ntypes; i++) { + if (!atom->radius_consistency(i,radtype)) + error->all(FLERR,"Pair lubricateU requires monodisperse particles"); + if (i > 1 && radtype != rad) + error->all(FLERR,"Pair lubricateU requires monodisperse particles"); + rad = radtype; + } + + // set the isotropic constants depending on the volume fraction + // vol_T = total volume + + double vol_T = domain->xprd*domain->yprd*domain->zprd; + + // assuming monodisperse spheres, vol_P = volume of the particles + + double tmp = 0.0; + if (atom->radius) tmp = atom->radius[0]; + double radi; + MPI_Allreduce(&tmp,&radi,1,MPI_DOUBLE,MPI_MAX,world); + + double vol_P = atom->natoms * (4.0/3.0)*MY_PI*pow(radi,3); + + // vol_f = volume fraction + + double vol_f = vol_P/vol_T; + + // set the isotropic constant + + if (flaglog == 0) { + R0 = 6*MY_PI*mu*radi*(1.0 + 2.16*vol_f); + RT0 = 8*MY_PI*mu*pow(radi,3); // not actually needed + RS0 = 20.0/3.0*MY_PI*mu*pow(radi,3)*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f); + } else { + R0 = 6*MY_PI*mu*radi*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f); + RT0 = 8*MY_PI*mu*pow(radi,3)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); + RS0 = 20.0/3.0*MY_PI*mu*pow(radi,3)*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f); + } +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairLubricateU::init_one(int i, int j) +{ + if (setflag[i][j] == 0) { + cut_inner[i][j] = mix_distance(cut_inner[i][i],cut_inner[j][j]); + cut[i][j] = mix_distance(cut[i][i],cut[j][j]); + } + + cut_inner[j][i] = cut_inner[i][j]; + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLubricateU::write_restart(FILE *fp) +{ + write_restart_settings(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(&cut_inner[i][j],sizeof(double),1,fp); + fwrite(&cut[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLubricateU::read_restart(FILE *fp) +{ + read_restart_settings(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(&cut_inner[i][j],sizeof(double),1,fp); + fread(&cut[i][j],sizeof(double),1,fp); + } + MPI_Bcast(&cut_inner[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLubricateU::write_restart_settings(FILE *fp) +{ + fwrite(&mu,sizeof(double),1,fp); + fwrite(&flaglog,sizeof(int),1,fp); + fwrite(&cut_inner_global,sizeof(double),1,fp); + fwrite(&cut_global,sizeof(double),1,fp); + fwrite(&offset_flag,sizeof(int),1,fp); + fwrite(&mix_flag,sizeof(int),1,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLubricateU::read_restart_settings(FILE *fp) +{ + int me = comm->me; + if (me == 0) { + fread(&mu,sizeof(double),1,fp); + fread(&flaglog,sizeof(int),1,fp); + fread(&cut_inner_global,sizeof(double),1,fp); + fread(&cut_global,sizeof(double),1,fp); + fread(&offset_flag,sizeof(int),1,fp); + fread(&mix_flag,sizeof(int),1,fp); + } + MPI_Bcast(&mu,1,MPI_DOUBLE,0,world); + MPI_Bcast(&flaglog,1,MPI_INT,0,world); + MPI_Bcast(&cut_inner_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&offset_flag,1,MPI_INT,0,world); + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); +} + +/*---------------------------------------------------------------------------*/ + +void PairLubricateU::copy_vec_uo(int inum, double *xcg, + double **v, double **omega) +{ + int i,j,ii; + int *ilist; + int itype; + double radi; + double inertia; + + double *rmass = atom->rmass; + int *type = atom->type; + + ilist = list->ilist; + + for (ii=0;iiradius[i]; + inertia = 0.4*rmass[i]*radi*radi; + + for (j=0;j<3;j++) { + v[i][j] = xcg[6*ii+j]; + omega[i][j] = xcg[6*ii+j+3]; + } + } +} + +/*---------------------------------------------------------------------------*/ + +void PairLubricateU::copy_uo_vec(int inum, double **f, double **torque, + double *RU) +{ + int i,j,ii; + int *ilist; + + ilist = list->ilist; + + for (ii=0;iiv; + double **omega = atom->omega; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + buf[m++] = omega[j][0]; + buf[m++] = omega[j][1]; + buf[m++] = omega[j][2]; + } + return 6; +} + +/* ---------------------------------------------------------------------- */ + +void PairLubricateU::unpack_comm(int n, int first, double *buf) +{ + int i,m,last; + + double **v = atom->v; + double **omega = atom->omega; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + omega[i][0] = buf[m++]; + omega[i][1] = buf[m++]; + omega[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + +double PairLubricateU::dot_vec_vec(int N, double *x, double *y) +{ + double dotp=0.0; + for (int i = 0; i < N; i++) dotp += x[i]*y[i]; + return dotp; +} diff --git a/src/FLD/pair_lubricateU.h b/src/FLD/pair_lubricateU.h new file mode 100644 index 0000000000..49bd43b9db --- /dev/null +++ b/src/FLD/pair_lubricateU.h @@ -0,0 +1,74 @@ +/* ---------------------------------------------------------------------- + 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(lubricateU,PairLubricateU) + +#else + +#ifndef LMP_PAIR_LUBRICATEU_H +#define LMP_PAIR_LUBRICATEU_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairLubricateU : public Pair { + public: + PairLubricateU(class LAMMPS *); + virtual ~PairLubricateU(); + virtual void compute(int, int); + virtual void settings(int, char **); + void coeff(int, char **); + double init_one(int, int); + virtual void init_style(); + void write_restart(FILE *); + void read_restart(FILE *); + void write_restart_settings(FILE *); + void read_restart_settings(FILE *); + int pack_comm(int, int *, double *, int, int *); + void unpack_comm(int, int, double *); + + protected: + double cut_inner_global,cut_global; + double mu; + int flaglog; + double gdot,Ef[3][3]; + double **cut_inner,**cut; + void allocate(); + double R0,RT0,RS0; + + int nmax; + double **fl,**Tl,**xl; + + int cgmax; + double *bcg,*xcg,*rcg,*rcg1,*pcg,*RU; + + void compute_RE(); + virtual void compute_RE(double **); + void compute_RU(); + virtual void compute_RU(double **); + virtual void compute_Fh(double **); + void stage_one(); + void intermediates(int, double **); + void stage_two(double **); + void copy_vec_uo(int, double *, double **, double **); + void copy_uo_vec(int, double **, double **, double *); + double dot_vec_vec(int , double *, double *); +}; + +} + +#endif +#endif diff --git a/src/FLD/pair_lubricateU_poly.cpp b/src/FLD/pair_lubricateU_poly.cpp new file mode 100644 index 0000000000..2fbdd0c730 --- /dev/null +++ b/src/FLD/pair_lubricateU_poly.cpp @@ -0,0 +1,1134 @@ +/* ---------------------------------------------------------------------- + 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: Amit Kumar and Michael Bybee (UIUC) + Pieter in 't Veld (BASF), code restructuring + Dave Heine (Corning), polydispersity +------------------------------------------------------------------------- */ + +#include "mpi.h" +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_lubricateU_poly.h" +#include "atom.h" +#include "atom_vec.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "domain.h" +#include "update.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +#define TOL 1E-3 // tolerance for conjugate gradient + +/* ---------------------------------------------------------------------- */ + +PairLubricateUPoly::PairLubricateUPoly(LAMMPS *lmp) : + PairLubricateU(lmp) {} + +/* ---------------------------------------------------------------------- + It first has to solve for the velocity of the particles such that + the net force on the particles is zero. NOTE: it has to be the last + type of pair interaction specified in the input file. Also, it + assumes that no other types of interactions, like k-space, is + present. As already mentioned, the net force on the particles after + this pair interaction would be identically zero. + ---------------------------------------------------------------------- */ + +void PairLubricateUPoly::compute(int eflag, int vflag) +{ + int i,j; + + int nlocal = atom->nlocal; + int nghost = atom->nghost; + int nall = nlocal + nghost; + + double **x = atom->x; + double **f = atom->f; + double **torque = atom->torque; + + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + // grow per-atom arrays if necessary + // need to be atom->nmax in length + + if (atom->nmax > nmax) { + memory->destroy(fl); + memory->destroy(Tl); + memory->destroy(xl); + nmax = atom->nmax; + memory->create(fl,nmax,3,"pair:fl"); + memory->create(Tl,nmax,3,"pair:Tl"); + memory->create(xl,nmax,3,"pair:xl"); + } + + if (6*list->inum > cgmax) { + memory->sfree(bcg); + memory->sfree(xcg); + memory->sfree(rcg); + memory->sfree(rcg1); + memory->sfree(pcg); + memory->sfree(RU); + cgmax = 6*list->inum; + memory->create(bcg,cgmax,"pair:bcg"); + memory->create(xcg,cgmax,"pair:bcg"); + memory->create(rcg,cgmax,"pair:bcg"); + memory->create(rcg1,cgmax,"pair:bcg"); + memory->create(pcg,cgmax,"pair:bcg"); + memory->create(RU,cgmax,"pair:bcg"); + } + + // Added to implement midpoint integration scheme + // Save force, torque found so far. Also save the positions + + for (i=0;ix,1); + + // Find positions at half the timestep and store in xl + + intermediates(nall,xl); + + // Store back the saved forces and torques in original arrays + + for(i=0;iinum; + int *ilist = list->ilist; + int *type = atom->type; + int newton_pair = force->newton_pair; + + double alpha,beta; + double normi,error,normig; + double send[2],recv[2],rcg_dot_rcg; + double mo_inertia,radi; + + double **v = atom->v; + double **f = atom->f; + double **omega = atom->omega; + double **angmom = atom->angmom; + double **torque = atom->torque; + double *radius = atom->radius; + + // First compute R_FE*E + + compute_RE(x); + + // Reverse communication of forces and torques to + // accumulate the net force on each of the particles + + if (newton_pair) comm->reverse_comm(); + + // CONJUGATE GRADIENT + // Find the right hand side= -ve of all forces/torques + // b = 6*Npart in overall size + + for(ii = 0; ii < inum; ii++) { + i = ilist[ii]; + for (j = 0; j < 3; j++) { + bcg[6*ii+j] = -f[i][j]; + bcg[6*ii+j+3] = -torque[i][j]; + } + } + + // Start solving the equation : F^H = -F^P -F^B - F^H_{Ef} + // Store initial guess for velocity and angular-velocities/angular momentum + // NOTE velocities and angular velocities are assumed relative to the fluid + + for (ii=0;iiforward_comm_pair(this); + + // Find initial residual + + compute_RU(x); + + // reverse communication of forces and torques + + if (newton_pair) comm->reverse_comm(); + + copy_uo_vec(inum,f,torque,RU); + + for (i=0;i<6*inum;i++) + rcg[i] = bcg[i] - RU[i]; + + // Set initial conjugate direction + + for (i=0;i<6*inum;i++) + pcg[i] = rcg[i]; + + // Find initial norm of the residual or norm of the RHS (either is fine) + + normi = dot_vec_vec(6*inum,bcg,bcg); + + MPI_Allreduce(&normi,&normig,1,MPI_DOUBLE,MPI_SUM,world); + + // Loop until convergence + + do { + // find R*p + + copy_vec_uo(inum,pcg,v,omega); + + // set velocities for ghost particles + + comm->forward_comm_pair(this); + + compute_RU(x); + + // reverse communication of forces and torques + + if (newton_pair) comm->reverse_comm(); + + copy_uo_vec(inum,f,torque,RU); + + // Find alpha + + send[0] = dot_vec_vec(6*inum,rcg,rcg); + send[1] = dot_vec_vec(6*inum,RU,pcg); + + MPI_Allreduce(send,recv,2,MPI_DOUBLE,MPI_SUM,world); + + alpha = recv[0]/recv[1]; + rcg_dot_rcg = recv[0]; + + // Find new x + + for (i=0;i<6*inum;i++) + xcg[i] = xcg[i] + alpha*pcg[i]; + + // find new residual + + for (i=0;i<6*inum;i++) + rcg1[i] = rcg[i] - alpha*RU[i]; + + // find beta + + send[0] = dot_vec_vec(6*inum,rcg1,rcg1); + + MPI_Allreduce(send,recv,1,MPI_DOUBLE,MPI_SUM,world); + + beta = recv[0]/rcg_dot_rcg; + + // Find new conjugate direction + + for (i=0;i<6*inum;i++) + pcg[i] = rcg1[i] + beta*pcg[i]; + + for (i=0;i<6*inum;i++) + rcg[i] = rcg1[i]; + + // Find relative error + + error = sqrt(recv[0]/normig); + + } while (error > TOL); + + + // update the final converged velocities in respective arrays + + copy_vec_uo(inum,xcg,v,omega); + + // set velocities for ghost particles + + comm->forward_comm_pair(this); + + // compute the viscosity/pressure + + if (evflag && stage == 2) compute_Fh(x); + + // find actual particle's velocities from relative velocities + // only non-zero component of fluid's vel : vx=gdot*y and wz=-gdot/2 + + for (ii=0;iitype; + int nlocal = atom->nlocal; + int nghost = atom->nghost; + int newton_pair = force->newton_pair; + int overlaps = 0; + + double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz; + double rsq,r,h_sep,radi,radj; + double vr1,vr2,vr3,vnnr,vn1,vn2,vn3; + double vt1,vt2,vt3; + double inv_inertia; + double vi[3],vj[3],wi[3],wj[3],xl[3],jl[3],pre[2]; + + double **v = atom->v; + double **f = atom->f; + double **omega = atom->omega; + double **angmom = atom->angmom; + double **torque = atom->torque; + double *radius = atom->radius; + + double beta[2][5]; + double vxmu2f = force->vxmu2f; + double a_sq = 0.0; + double a_sh = 0.0; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + beta[0][0] = beta[1][0] = beta[1][4] = 0.0; + + // Set force to zero which is the final value after this pair interaction + + for (i=0;ireverse_comm(); // not really needed + + // Find additional contribution from the stresslets + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + radi = radius[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + pre[1] = 8.0*(pre[0] = MY_PI*mu*radi)*radi*radi; + pre[0] *= 6.0; + + // Find the contribution to stress from isotropic RS0 + // Set psuedo force to obtain the required contribution + // need to set delx and fy only + + fx = 0.0; delx = radi; + fy = vxmu2f*RS0*pow(radi,3)*gdot/2.0/radi; dely = 0.0; + fz = 0.0; delz = 0.0; + if (evflag) + ev_tally_xyz(i,i,nlocal,newton_pair,0.0,0.0,-fx,-fy,-fz,delx,dely,delz); + + // Find angular velocity + + wi[0] = omega[i][0]; + wi[1] = omega[i][1]; + wi[2] = omega[i][2]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + radj = radius[j]; + + if (rsq < cutsq[itype][jtype]) { + r = sqrt(rsq); + + // Use omega directly if it exists, else angmom + // angular momentum = I*omega = 2/5 * M*R^2 * omega + + wj[0] = omega[j][0]; + wj[1] = omega[j][1]; + wj[2] = omega[j][2]; + + // loc of the point of closest approach on particle i from its cente + // POC for j is in opposite direction as for i + + xl[0] = -delx/r*radi; + xl[1] = -dely/r*radi; + xl[2] = -delz/r*radi; + jl[0] = delx/r*radj; + jl[1] = dely/r*radj; + jl[2] = delz/r*radj; + + h_sep = r - radi-radj; + + // velocity at the point of closest approach on both particles + // v = v + omega_cross_xl + + // particle i + + vi[0] = v[i][0] + (wi[1]*xl[2] - wi[2]*xl[1]); + vi[1] = v[i][1] + (wi[2]*xl[0] - wi[0]*xl[2]); + vi[2] = v[i][2] + (wi[0]*xl[1] - wi[1]*xl[0]); + + // particle j + + vj[0] = v[j][0] + (wj[1]*jl[2] - wj[2]*jl[1]); + vj[1] = v[j][1] + (wj[2]*jl[0] - wj[0]*jl[2]); + vj[2] = v[j][2] + (wj[0]*jl[1] - wj[1]*jl[0]); + + + // Relative velocity at the point of closest approach + // include contribution from Einf of the fluid + + vr1 = vi[0] - vj[0] - + 2.0*(Ef[0][0]*xl[0] + Ef[0][1]*xl[1] + Ef[0][2]*xl[2]); + vr2 = vi[1] - vj[1] - + 2.0*(Ef[1][0]*xl[0] + Ef[1][1]*xl[1] + Ef[1][2]*xl[2]); + vr3 = vi[2] - vj[2] - + 2.0*(Ef[2][0]*xl[0] + Ef[2][1]*xl[1] + Ef[2][2]*xl[2]); + + // Normal component (vr.n)n + + vnnr = (vr1*delx + vr2*dely + vr3*delz)/r; + vn1 = vnnr*delx/r; + vn2 = vnnr*dely/r; + vn3 = vnnr*delz/r; + + // Tangential component vr - (vr.n)n + + vt1 = vr1 - vn1; + vt2 = vr2 - vn2; + vt3 = vr3 - vn3; + + // Find the scalar resistances a_sq, a_sh and a_pu + + // check for overlaps + + if (h_sep < 0.0) overlaps++; + + // If less than the minimum gap use the minimum gap instead + + if (r < cut_inner[itype][jtype]) + h_sep = cut_inner[itype][jtype] - radi-radj; + + // Scale h_sep by radi + + h_sep = h_sep/radi; + beta[0][1] = radj/radi; + beta[1][1] = 1.0 + beta[0][1]; + + /*beta0 = radj/radi; + beta1 = 1.0 + beta0;*/ + + // Scalar resistances + + if (flaglog) { + beta[0][2] = beta[0][1]*beta[0][1]; + beta[0][3] = beta[0][2]*beta[0][1]; + beta[0][4] = beta[0][3]*beta[0][1]; + beta[1][2] = beta[1][1]*beta[1][1]; + beta[1][3] = beta[1][2]*beta[1][1]; + double log_h_sep_beta13 = log(1.0/h_sep)/beta[1][3]; + double h_sep_beta11 = h_sep/beta[1][1]; + + a_sq = pre[0]*(beta[0][2]/beta[1][2]/h_sep + +((0.2+1.4*beta[0][1]+0.2*beta[0][2]) + +(1.0+18.0*(beta[0][1]+beta[0][3])-29.0*beta[0][2] + +beta[0][4])*h_sep_beta11/21.0)*log_h_sep_beta13); + + a_sh = pre[0]*((8.0*(beta[0][1]+beta[0][3])+4.0*beta[0][2])/15.0 + +(64.0-180.0*(beta[0][1]+beta[0][3])+232.0*beta[0][2] + +64.0*beta[0][4])*h_sep_beta11/375.0)*log_h_sep_beta13; + + /*a_sq = beta0*beta0/beta1/beta1/h_sep + +(1.0+7.0*beta0+beta0*beta0)/5.0/pow(beta1,3)*log(1.0/h_sep); + a_sq += (1.0+18.0*beta0-29.0*beta0*beta0+18.0*pow(beta0,3) + +pow(beta0,4))/21.0/pow(beta1,4)*h_sep*log(1.0/h_sep); + a_sq *= 6.0*MY_PI*mu*radi; + + a_sh = 4.0*beta0*(2.0+beta0 + +2.0*beta0*beta0)/15.0/pow(beta1,3)*log(1.0/h_sep); + a_sh += 4.0*(16.0-45.0*beta0+58.0*beta0*beta0-45.0*pow(beta0,3) + +16.0*pow(beta0,4))/375.0/pow(beta1,4)*h_sep*log(1.0/h_sep); + a_sh *= 6.0*MY_PI*mu*radi;*/ + } else { + //a_sq = 6.0*MY_PI*mu*radi*(beta0*beta0/beta1/beta1/h_sep); + a_sq = pre[0]*(beta[0][1]*beta[0][1]/(beta[1][1]*beta[1][1]*h_sep)); + } + + // Find force due to squeeze type motion + + fx = a_sq*vn1; + fy = a_sq*vn2; + fz = a_sq*vn3; + + // Find force due to all shear kind of motions + + if (flaglog) { + fx = fx + a_sh*vt1; + fy = fy + a_sh*vt2; + fz = fz + a_sh*vt3; + } + + // Scale forces to obtain in appropriate units + + fx = vxmu2f*fx; + fy = vxmu2f*fy; + fz = vxmu2f*fz; + + // set j = nlocal so that only I gets tallied + + if (evflag) ev_tally_xyz(i,nlocal,nlocal,0, + 0.0,0.0,-fx,-fy,-fz,delx,dely,delz); + } + } + } +} + +/* ---------------------------------------------------------------------- + computes R_FU * U +---------------------------------------------------------------------- */ + +void PairLubricateUPoly::compute_RU(double **x) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + int *ilist,*jlist,*numneigh,**firstneigh; + + int *type = atom->type; + int nlocal = atom->nlocal; + int nghost = atom->nghost; + int overlaps = 0; + + double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,tx,ty,tz; + double rsq,r,radi,radj,h_sep; + //double beta0,beta1; + double vr1,vr2,vr3,vnnr,vn1,vn2,vn3; + double vt1,vt2,vt3,wdotn,wt1,wt2,wt3; + double inv_inertia; + double vi[3],vj[3],wi[3],wj[3],xl[3],jl[3],pre[2]; + + double **v = atom->v; + double **f = atom->f; + double **omega = atom->omega; + double **angmom = atom->angmom; + double **torque = atom->torque; + double *radius = atom->radius; + + double beta[2][5]; + double vxmu2f = force->vxmu2f; + double a_sq = 0.0; + double a_sh = 0.0; + double a_pu = 0.0; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + beta[0][0] = beta[1][0] = beta[1][4] = 0.0; + + // Initialize f to zero + + for (i=0;itype; + int overlaps = 0; + + double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,tx,ty,tz; + double rsq,r,h_sep,radi,radj; + //double beta0,beta1,lhsep; + double vr1,vr2,vr3,vnnr,vn1,vn2,vn3; + double vt1,vt2,vt3; + double xl[3],pre[2]; + + double **f = atom->f; + double **torque = atom->torque; + double *radius = atom->radius; + + double beta[2][5]; + double vxmu2f = force->vxmu2f; + double a_sq = 0.0; + double a_sh = 0.0; + double a_pu = 0.0; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + beta[0][0] = beta[1][0] = beta[1][4] = 0.0; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + radi = radius[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + pre[1] = 8.0*(pre[0] = MY_PI*mu*radi)*radi*radi; + pre[0] *= 6.0; + + // No contribution from isotropic terms due to E + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + radj = radius[j]; + + if (rsq < cutsq[itype][jtype]) { + r = sqrt(rsq); + + // loc of the point of closest approach on particle i from its center + + xl[0] = -delx/r*radi; + xl[1] = -dely/r*radi; + xl[2] = -delz/r*radi; + + // Find the scalar resistances a_sq and a_sh + + h_sep = r - radi-radj; + + // check for overlaps + + if (h_sep < 0.0) overlaps++; + + // If less than the minimum gap use the minimum gap instead + + if (r < cut_inner[itype][jtype]) + h_sep = cut_inner[itype][jtype] - radi-radj; + + // Scale h_sep by radi + + h_sep = h_sep/radi; + beta[0][1] = radj/radi; + beta[1][1] = 1.0 + beta[0][1]; + + /*beta0 = radj/radi; + beta1 = 1.0 + beta0; + lhsep = log(1.0/h_sep);*/ + + // Scalar resistance for Squeeze type motions + + + if (flaglog) { + beta[0][2] = beta[0][1]*beta[0][1]; + beta[0][3] = beta[0][2]*beta[0][1]; + beta[0][4] = beta[0][3]*beta[0][1]; + beta[1][2] = beta[1][1]*beta[1][1]; + beta[1][3] = beta[1][2]*beta[1][1]; + double log_h_sep_beta13 = log(1.0/h_sep)/beta[1][3]; + double h_sep_beta11 = h_sep/beta[1][1]; + + a_sq = pre[0]*(beta[0][2]/beta[1][2]/h_sep + +((0.2+1.4*beta[0][1]+0.2*beta[0][2]) + +(1.0+18.0*(beta[0][1]+beta[0][3])-29.0*beta[0][2] + +beta[0][4])*h_sep_beta11/21.0)*log_h_sep_beta13); + + a_sh = pre[0]*((8.0*(beta[0][1]+beta[0][3])+4.0*beta[0][2])/15.0 + +(64.0-180.0*(beta[0][1]+beta[0][3])+232.0*beta[0][2] + +64.0*beta[0][4])*h_sep_beta11/375.0)*log_h_sep_beta13; + + a_pu = pre[1]*((0.4*beta[0][1]+0.1*beta[0][2])*beta[1][1] + +(0.128-0.132*beta[0][1]+0.332*beta[0][2] + +0.172*beta[0][3])*h_sep)*log_h_sep_beta13; + + /*//a_sq = 6*MY_PI*mu*radi*(1.0/4.0/h_sep + 9.0/40.0*log(1/h_sep)); + a_sq = beta0*beta0/beta1/beta1/h_sep + +(1.0+7.0*beta0+beta0*beta0)/5.0/pow(beta1,3)*lhsep; + a_sq += (1.0+18.0*beta0-29.0*beta0*beta0+18.0*pow(beta0,3) + +pow(beta0,4))/21.0/pow(beta1,4)*h_sep*lhsep; + a_sq *= 6.0*MY_PI*mu*radi; + + a_sh = 4.0*beta0*(2.0+beta0 + +2.0*beta0*beta0)/15.0/pow(beta1,3)*log(1.0/h_sep); + a_sh += 4.0*(16.0-45.0*beta0+58.0*beta0*beta0-45.0*pow(beta0,3) + +16.0*pow(beta0,4))/375.0/pow(beta1,4)*h_sep*log(1.0/h_sep); + a_sh *= 6.0*MY_PI*mu*radi; + + a_pu = beta0*(4.0+beta0)/10.0/beta1/beta1*log(1.0/h_sep); + a_pu += (32.0-33.0*beta0+83.0*beta0*beta0 + +43.0*pow(beta0,3))/250.0/pow(beta1,3)*h_sep*log(1.0/h_sep); + a_pu *= 8.0*MY_PI*mu*pow(radi,3);*/ + } else + a_sq = pre[0]*(beta[0][1]*beta[0][1]/(beta[1][1]*beta[1][1]*h_sep)); + + // Relative velocity at the point of closest approach due to Ef only + + vr1 = -2.0*(Ef[0][0]*xl[0] + Ef[0][1]*xl[1] + Ef[0][2]*xl[2]); + vr2 = -2.0*(Ef[1][0]*xl[0] + Ef[1][1]*xl[1] + Ef[1][2]*xl[2]); + vr3 = -2.0*(Ef[2][0]*xl[0] + Ef[2][1]*xl[1] + Ef[2][2]*xl[2]); + + // Normal component (vr.n)n + + vnnr = (vr1*delx + vr2*dely + vr3*delz)/r; + vn1 = vnnr*delx/r; + vn2 = vnnr*dely/r; + vn3 = vnnr*delz/r; + + // Tangential component vr - (vr.n)n + + vt1 = vr1 - vn1; + vt2 = vr2 - vn2; + vt3 = vr3 - vn3; + + // Find force due to squeeze type motion + + fx = a_sq*vn1; + fy = a_sq*vn2; + fz = a_sq*vn3; + + // Find force due to all shear kind of motions + + if (flaglog) { + fx = fx + a_sh*vt1; + fy = fy + a_sh*vt2; + fz = fz + a_sh*vt3; + } + + // Scale forces to obtain in appropriate units + + fx = vxmu2f*fx; + fy = vxmu2f*fy; + fz = vxmu2f*fz; + + // Add to the total forc + + f[i][0] -= fx; + f[i][1] -= fy; + f[i][2] -= fz; + + // Find torque due to this force + + if (flaglog) { + tx = xl[1]*fz - xl[2]*fy; + ty = xl[2]*fx - xl[0]*fz; + tz = xl[0]*fy - xl[1]*fx; + + // Why a scale factor ? + + torque[i][0] -= vxmu2f*tx; + torque[i][1] -= vxmu2f*ty; + torque[i][2] -= vxmu2f*tz; + + // NOTE No a_pu term needed as they add up to zero + } + } + } + } + + int print_overlaps = 0; + if (print_overlaps && overlaps) + printf("Number of overlaps=%d\n",overlaps); +} + +/*----------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairLubricateUPoly::settings(int narg, char **arg) +{ + double vol_T; + int itype; + + if (narg != 5) error->all(FLERR,"Illegal pair_style command"); + + mu = atof(arg[0]); + flaglog = atoi(arg[1]); + cut_inner_global = atof(arg[2]); + cut_global = atof(arg[3]); + gdot = atof(arg[4]); + + // reset cutoffs that have been explicitly set + + if (allocated) { + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i+1; j <= atom->ntypes; j++) + if (setflag[i][j]) { + cut_inner[i][j] = cut_inner_global; + cut[i][j] = cut_global; + } + } + + // Store the rate of strain tensor + + Ef[0][0] = 0.0; + Ef[0][1] = 0.5*gdot; + Ef[0][2] = 0.0; + Ef[1][0] = 0.5*gdot; + Ef[1][1] = 0.0; + Ef[1][2] = 0.0; + Ef[2][0] = 0.0; + Ef[2][1] = 0.0; + Ef[2][2] = 0.0; + + // Set the isotropic constants depending on the volume fraction + + // Find the total volume + + vol_T = domain->xprd*domain->yprd*domain->zprd; + + // Assuming monodisperse spheres, find the volume of the particles + + int nlocal = atom->nlocal; + int *type = atom->type; + + double volP = 0.0; + for (int i = 0; i < nlocal; i++) + volP += (4.0/3.0)*MY_PI*pow(atom->radius[i],3); + double vol_P; + MPI_Allreduce(&volP,&vol_P,1,MPI_DOUBLE,MPI_SUM,world); + double vol_f = vol_P/vol_T; + + //DRH volume fraction needs to be defined manually + // if excluded volume regions are present + + vol_f=0.5; + if (!comm->me) { + if(logfile) + fprintf(logfile, "lubricateU: vol_f = %g, vol_p = %g, vol_T = %g\n", + vol_f,vol_P,vol_T); + if (screen) + fprintf(screen, "lubricateU: vol_f = %g, vol_p = %g, vol_T = %g\n", + vol_f,vol_P,vol_T); + } + + // Set the isotropic constant + + if (flaglog == 0) { + R0 = 6*MY_PI*mu*(1.0 + 2.16*vol_f); + RT0 = 8*MY_PI*mu; // Not needed actually + RS0 = 20.0/3.0*MY_PI*mu*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f); + } else { + R0 = 6*MY_PI*mu*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f); + RT0 = 8*MY_PI*mu*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); + RS0 = 20.0/3.0*MY_PI*mu*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f); + } +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairLubricateUPoly::init_style() +{ + if (force->newton_pair == 1) + error->all(FLERR,"Pair lubricateU/poly requires newton pair off"); + if (comm->ghost_velocity == 0) + error->all(FLERR, + "Pair lubricateU/poly requires ghost atoms store velocity"); + if (!atom->sphere_flag) + error->all(FLERR,"Pair lubricate/poly requires atom style sphere"); + + // insure all particles are finite-size + // for pair hybrid, should limit test to types using the pair style + + double *radius = atom->radius; + int *type = atom->type; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (radius[i] == 0.0) + error->one(FLERR,"Pair lubricate/poly requires extended particles"); + + int irequest = neighbor->request(this); + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->full = 1; +} diff --git a/src/FLD/pair_lubricateU_poly.h b/src/FLD/pair_lubricateU_poly.h new file mode 100644 index 0000000000..4596a13da5 --- /dev/null +++ b/src/FLD/pair_lubricateU_poly.h @@ -0,0 +1,45 @@ +/* ---------------------------------------------------------------------- + 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(lubricateU/poly,PairLubricateUPoly) + +#else + +#ifndef LMP_PAIR_LUBRICATEU_POLY_H +#define LMP_PAIR_LUBRICATEU_POLY_H + +#include "pair_lubricateU.h" + +namespace LAMMPS_NS { + +class PairLubricateUPoly : public PairLubricateU { + public: + PairLubricateUPoly(class LAMMPS *); + ~PairLubricateUPoly() {} + void compute(int, int); + void settings(int, char **); + void init_style(); + + private: + void iterate(double **, int); + void compute_RE(double **); + void compute_RU(double **); + void compute_Fh(double **); +}; + +} + +#endif +#endif diff --git a/src/FLD/pair_lubricate_poly.cpp b/src/FLD/pair_lubricate_poly.cpp new file mode 100644 index 0000000000..1fa5d885ef --- /dev/null +++ b/src/FLD/pair_lubricate_poly.cpp @@ -0,0 +1,459 @@ +/* ---------------------------------------------------------------------- + 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: Randy Schunk (SNL) + Amit Kumar and Michael Bybee (UIUC) + Dave Heine (Corning), polydispersity +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_lubricate_poly.h" +#include "atom.h" +#include "atom_vec.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "domain.h" +#include "update.h" +#include "modify.h" +#include "fix.h" +#include "fix_deform.h" +#include "memory.h" +#include "random_mars.h" +#include "math_const.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +// same as fix_deform.cpp + +enum{NO_REMAP,X_REMAP,V_REMAP}; + +/* ---------------------------------------------------------------------- */ + +PairLubricatePoly::PairLubricatePoly(LAMMPS *lmp) : PairLubricate(lmp) +{ + no_virial_fdotr_compute = 1; +} + +/* ---------------------------------------------------------------------- */ + +void PairLubricatePoly::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,fpair,fx,fy,fz,tx,ty,tz; + double rsq,r,h_sep,h_sepj,beta0,beta1,betaj,betaj1,radi,radj,tfmag; + double vr1,vr2,vr3,vnnr,vn1,vn2,vn3; + double vt1,vt2,vt3,wt1,wt2,wt3,wdotn; + double inertia,inv_inertia,vRS0; + double vi[3],vj[3],wi[3],wj[3],xl[3],jl[3]; + double a_sq,a_sh,a_pu,Fbmag,del,delmin,eta; + int *ilist,*jlist,*numneigh,**firstneigh; + double lamda[3],vstream[3]; + + double vxmu2f = force->vxmu2f; + + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + double **omega = atom->omega; + double **angmom = atom->angmom; + double **torque = atom->torque; + double *radius = atom->radius; + double *mass = atom->mass; + double *rmass = atom->rmass; + int *type = atom->type; + int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; + + int overlaps = 0; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // subtract streaming component of velocity, omega, angmom + // assume fluid streaming velocity = box deformation rate + // vstream = (ux,uy,uz) + // ux = h_rate[0]*x + h_rate[5]*y + h_rate[4]*z + // uy = h_rate[1]*y + h_rate[3]*z + // uz = h_rate[2]*z + // omega_new = omega - curl(vstream)/2 + // angmom_new = angmom - I*curl(vstream)/2 + // Ef = (grad(vstream) + (grad(vstream))^T) / 2 + + if (shearing) { + double *h_rate = domain->h_rate; + double *h_ratelo = domain->h_ratelo; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + itype = type[i]; + radi = radius[i]; + domain->x2lamda(x[i],lamda); + vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + + h_rate[4]*lamda[2] + h_ratelo[0]; + vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; + vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2]; + v[i][0] -= vstream[0]; + v[i][1] -= vstream[1]; + v[i][2] -= vstream[2]; + + omega[i][0] += 0.5*h_rate[3]; + omega[i][1] -= 0.5*h_rate[4]; + omega[i][2] += 0.5*h_rate[5]; + } + + // set Ef from h_rate in strain units + + Ef[0][0] = h_rate[0]/domain->xprd; + Ef[1][1] = h_rate[1]/domain->yprd; + Ef[2][2] = h_rate[2]/domain->zprd; + Ef[0][1] = Ef[1][0] = 0.5 * h_rate[5]/domain->yprd; + Ef[0][2] = Ef[2][0] = 0.5 * h_rate[4]/domain->zprd; + Ef[1][2] = Ef[2][1] = 0.5 * h_rate[3]/domain->zprd; + + // copy updated omega to the ghost particles + // no need to do this if not shearing since comm->ghost_velocity is set + + comm->forward_comm_pair(this); + } + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + radi = radius[i]; + + // angular velocity + + wi[0] = omega[i][0]; + wi[1] = omega[i][1]; + wi[2] = omega[i][2]; + + // FLD contribution to force and torque due to isotropic terms + // FLD contribution to stress from isotropic RS0 + + if (flagfld) { + f[i][0] -= vxmu2f*R0*radi*v[i][0]; + f[i][1] -= vxmu2f*R0*radi*v[i][1]; + f[i][2] -= vxmu2f*R0*radi*v[i][2]; + torque[i][0] -= vxmu2f*RT0*pow(radi,3)*wi[0]; + torque[i][1] -= vxmu2f*RT0*pow(radi,3)*wi[1]; + torque[i][2] -= vxmu2f*RT0*pow(radi,3)*wi[2]; + + if (shearing && vflag_either) { + vRS0 = -vxmu2f * RS0*pow(radi,3); + v_tally_tensor(i,i,nlocal,newton_pair, + vRS0*Ef[0][0],vRS0*Ef[1][1],vRS0*Ef[2][2], + vRS0*Ef[0][1],vRS0*Ef[0][2],vRS0*Ef[1][2]); + } + } + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + radj = atom->radius[j]; + + if (rsq < cutsq[itype][jtype]) { + r = sqrt(rsq); + + // angular momentum = I*omega = 2/5 * M*R^2 * omega + + wj[0] = omega[j][0]; + wj[1] = omega[j][1]; + wj[2] = omega[j][2]; + + // xl = point of closest approach on particle i from its center + + xl[0] = -delx/r*radi; + xl[1] = -dely/r*radi; + xl[2] = -delz/r*radi; + jl[0] = -delx/r*radj; + jl[1] = -dely/r*radj; + jl[2] = -delz/r*radj; + + // velocity at the point of closest approach on both particles + // v = v + omega_cross_xl - Ef.xl + + // particle i + + vi[0] = v[i][0] + (wi[1]*xl[2] - wi[2]*xl[1]) + - (Ef[0][0]*xl[0] + Ef[0][1]*xl[1] + Ef[0][2]*xl[2]); + + vi[1] = v[i][1] + (wi[2]*xl[0] - wi[0]*xl[2]) + - (Ef[1][0]*xl[0] + Ef[1][1]*xl[1] + Ef[1][2]*xl[2]); + + vi[2] = v[i][2] + (wi[0]*xl[1] - wi[1]*xl[0]) + - (Ef[2][0]*xl[0] + Ef[2][1]*xl[1] + Ef[2][2]*xl[2]); + + // particle j + + vj[0] = v[j][0] - (wj[1]*jl[2] - wj[2]*jl[1]) + + (Ef[0][0]*jl[0] + Ef[0][1]*jl[1] + Ef[0][2]*jl[2]); + + vj[1] = v[j][1] - (wj[2]*jl[0] - wj[0]*jl[2]) + + (Ef[1][0]*jl[0] + Ef[1][1]*jl[1] + Ef[1][2]*jl[2]); + + vj[2] = v[j][2] - (wj[0]*jl[1] - wj[1]*jl[0]) + + (Ef[2][0]*jl[0] + Ef[2][1]*jl[1] + Ef[2][2]*jl[2]); + + // scalar resistances XA and YA + + h_sep = r - radi-radj; + + // check for overlaps + + if (h_sep < 0.0) overlaps++; + + // if less than the minimum gap use the minimum gap instead + + if (r < cut_inner[itype][jtype]) + h_sep = cut_inner[itype][jtype] - radi-radj; + + // scale h_sep by radi + + h_sep = h_sep/radi; + beta0 = radj/radi; + beta1 = 1.0 + beta0; + + // scalar resistances + + if (flaglog) { + a_sq = beta0*beta0/beta1/beta1/h_sep + + (1.0+7.0*beta0+beta0*beta0)/5.0/pow(beta1,3)*log(1.0/h_sep); + a_sq += (1.0+18.0*beta0-29.0*beta0*beta0+18.0 * + pow(beta0,3)+pow(beta0,4))/21.0/pow(beta1,4) * + h_sep*log(1.0/h_sep); + a_sq *= 6.0*MY_PI*mu*radi; + a_sh = 4.0*beta0*(2.0+beta0+2.0*beta0*beta0)/15.0/pow(beta1,3) * + log(1.0/h_sep); + a_sh += 4.0*(16.0-45.0*beta0+58.0*beta0*beta0-45.0*pow(beta0,3) + + 16.0*pow(beta0,4))/375.0/pow(beta1,4) * + h_sep*log(1.0/h_sep); + a_sh *= 6.0*MY_PI*mu*radi; + a_pu = beta0*(4.0+beta0)/10.0/beta1/beta1*log(1.0/h_sep); + a_pu += (32.0-33.0*beta0+83.0*beta0*beta0+43.0 * + pow(beta0,3))/250.0/pow(beta1,3)*h_sep*log(1.0/h_sep); + a_pu *= 8.0*MY_PI*mu*pow(radi,3); + } else a_sq = 6.0*MY_PI*mu*radi*(beta0*beta0/beta1/beta1/h_sep); + + // relative velocity at the point of closest approach + // includes fluid velocity + + vr1 = vi[0] - vj[0]; + vr2 = vi[1] - vj[1]; + vr3 = vi[2] - vj[2]; + + // normal component (vr.n)n + + vnnr = (vr1*delx + vr2*dely + vr3*delz)/r; + vn1 = vnnr*delx/r; + vn2 = vnnr*dely/r; + vn3 = vnnr*delz/r; + + // tangential component vr - (vr.n)n + + vt1 = vr1 - vn1; + vt2 = vr2 - vn2; + vt3 = vr3 - vn3; + + // force due to squeeze type motion + + fx = a_sq*vn1; + fy = a_sq*vn2; + fz = a_sq*vn3; + + // force due to all shear kind of motions + + if (flaglog) { + fx = fx + a_sh*vt1; + fy = fy + a_sh*vt2; + fz = fz + a_sh*vt3; + } + + // scale forces for appropriate units + + fx *= vxmu2f; + fy *= vxmu2f; + fz *= vxmu2f; + + // add to total force + + f[i][0] -= fx; + f[i][1] -= fy; + f[i][2] -= fz; + + // torque due to this force + + if (flaglog) { + tx = xl[1]*fz - xl[2]*fy; + ty = xl[2]*fx - xl[0]*fz; + tz = xl[0]*fy - xl[1]*fx; + + torque[i][0] -= vxmu2f*tx; + torque[i][1] -= vxmu2f*ty; + torque[i][2] -= vxmu2f*tz; + + // torque due to a_pu + + wdotn = ((wi[0]-wj[0])*delx + (wi[1]-wj[1])*dely + + (wi[2]-wj[2])*delz)/r; + wt1 = (wi[0]-wj[0]) - wdotn*delx/r; + wt2 = (wi[1]-wj[1]) - wdotn*dely/r; + wt3 = (wi[2]-wj[2]) - wdotn*delz/r; + + tx = a_pu*wt1; + ty = a_pu*wt2; + tz = a_pu*wt3; + + torque[i][0] -= vxmu2f*tx; + torque[i][1] -= vxmu2f*ty; + torque[i][2] -= vxmu2f*tz; + + } + + // set j = nlocal so that only I gets tallied + + if (evflag) ev_tally_xyz(i,nlocal,nlocal,0, + 0.0,0.0,-fx,-fy,-fz,delx,dely,delz); + } + } + } + + // restore streaming component of velocity, omega, angmom + + if (shearing) { + double *h_rate = domain->h_rate; + double *h_ratelo = domain->h_ratelo; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + itype = type[i]; + radi = atom->radius[i]; + + domain->x2lamda(x[i],lamda); + vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + + h_rate[4]*lamda[2] + h_ratelo[0]; + vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; + vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2]; + v[i][0] += vstream[0]; + v[i][1] += vstream[1]; + v[i][2] += vstream[2]; + + omega[i][0] -= 0.5*h_rate[3]; + omega[i][1] += 0.5*h_rate[4]; + omega[i][2] -= 0.5*h_rate[5]; + } + } + + // to DEBUG: set print_overlaps to 1 + + int print_overlaps = 0; + if (print_overlaps) { + int overlaps_all; + MPI_Allreduce(&overlaps,&overlaps_all,1,MPI_INT,MPI_SUM,world); + if (overlaps_all && comm->me == 0) + printf("Number of overlaps = %d\n",overlaps); + } +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairLubricatePoly::init_style() +{ + if (force->newton_pair == 1) + error->all(FLERR,"Pair lubricate/poly requires newton pair off"); + if (comm->ghost_velocity == 0) + error->all(FLERR, + "Pair lubricate/poly requires ghost atoms store velocity"); + if (!atom->sphere_flag) + error->all(FLERR,"Pair lubricate/poly requires atom style sphere"); + + // insure all particles are finite-size + // for pair hybrid, should limit test to types using the pair style + + double *radius = atom->radius; + int *type = atom->type; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (radius[i] == 0.0) + error->one(FLERR,"Pair lubricate/poly requires extended particles"); + + int irequest = neighbor->request(this); + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->full = 1; + + // set the isotropic constants that depend on the volume fraction + // vol_T = total volume + + double vol_T = domain->xprd*domain->yprd*domain->zprd; + + double volP = 0.0; + for (int i = 0; i < nlocal; i++) + volP += (4.0/3.0)*MY_PI*pow(atom->radius[i],3); + double vol_P; + MPI_Allreduce(&volP,&vol_P,1,MPI_DOUBLE,MPI_SUM,world); + double vol_f = vol_P/vol_T; + + // set isotropic constants + + if (flaglog == 0) { + R0 = 6*MY_PI*mu*(1.0 + 2.16*vol_f); + RT0 = 8*MY_PI*mu; + RS0 = 20.0/3.0*MY_PI*mu*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f); + } else { + R0 = 6*MY_PI*mu*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f); + RT0 = 8*MY_PI*mu*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f); + RS0 = 20.0/3.0*MY_PI*mu*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f); + } + + // check for fix deform, if exists it must use "remap v" + + shearing = 0; + for (int i = 0; i < modify->nfix; i++) + if (strcmp(modify->fix[i]->style,"deform") == 0) { + shearing = 1; + if (((FixDeform *) modify->fix[i])->remapflag != V_REMAP) + error->all(FLERR,"Using pair lubricate/poly with inconsistent " + "fix deform remap option"); + } + + // set Ef = 0 since used whether shearing or not + + Ef[0][0] = Ef[0][1] = Ef[0][2] = 0.0; + Ef[1][0] = Ef[1][1] = Ef[1][2] = 0.0; + Ef[2][0] = Ef[2][1] = Ef[2][2] = 0.0; +} diff --git a/src/FLD/pair_lubricate_poly.h b/src/FLD/pair_lubricate_poly.h new file mode 100644 index 0000000000..e13edd3313 --- /dev/null +++ b/src/FLD/pair_lubricate_poly.h @@ -0,0 +1,38 @@ +/* ---------------------------------------------------------------------- + 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(lubricate/poly,PairLubricatePoly) + +#else + +#ifndef LMP_PAIR_LUBRICATE_POLY_H +#define LMP_PAIR_LUBRICATE_POLY_H + +#include "pair_lubricate.h" + +namespace LAMMPS_NS { + +class PairLubricatePoly : public PairLubricate { + public: + PairLubricatePoly(class LAMMPS *); + ~PairLubricatePoly() {} + void compute(int, int); + void init_style(); +}; + +} + +#endif +#endif diff --git a/src/Makefile b/src/Makefile index c0a9df822a..71d2066e05 100755 --- a/src/Makefile +++ b/src/Makefile @@ -13,7 +13,7 @@ OBJ = $(SRC:.cpp=.o) # Package variables -PACKAGE = asphere class2 colloid dipole gpu granular \ +PACKAGE = asphere class2 colloid dipole fld gpu granular \ kspace manybody mc meam molecule opt peri poems reax replica \ shock srd xtc