Homec4science
Diffusion lammps (lammm-master) Browse

lammps/src/FLD/pair_lubricateU_poly.cpp
a362ece3e7e1lammm-master

Last Change
Raw File

pair_lubricateU_poly.cpp
View Options

/* ----------------------------------------------------------------------
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;i<nlocal+nghost;i++) {
for (j=0;j<3;j++) {
fl[i][j] = f[i][j];
Tl[i][j] = torque[i][j];
xl[i][j] = x[i][j];
}
}
// Stage one of Midpoint method
// Solve for velocities based on intial positions
iterate(atom->x,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;i<nlocal+nghost;i++) {
for(j=0;j<3;j++) {
f[i][j] = fl[i][j];
torque[i][j] = Tl[i][j];
}
}
// Stage two: This will give the final velocities
iterate(xl,2);
}
/* ------------------------------------------------------------------------
Stage one of midpoint method
------------------------------------------------------------------------- */
void PairLubricateUPoly::iterate(double **x, int stage)
{
int i,j,ii,itype;
int inum = list->inum;
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;ii<inum;ii++)
for (j=0;j<3;j++) {
xcg[6*ii+j] = 0.0;
xcg[6*ii+j+3] = 0.0;
}
// Copy initial guess to the global arrays to be acted upon by R_{FU}
// and returned by f and torque arrays
copy_vec_uo(inum,xcg,v,omega);
// set velocities for ghost particles
comm->forward_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;ii<inum;ii++) {
i = ilist[ii];
itype = type[i];
radi = radius[i];
v[i][0] = v[i][0] + gdot*x[i][1];
omega[i][2] = omega[i][2] - gdot/2.0;
}
}
/* ------------------------------------------------------------------------
This function computes the final hydrodynamic force once the
velocities have converged.
------------------------------------------------------------------------- */
void PairLubricateUPoly::compute_Fh(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 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;i<nlocal+nghost;i++)
for (j=0;j<3;j++) {
f[i][j] = 0.0;
torque[i][j] = 0.0;
}
// reverse communication of forces and torques
if (newton_pair) comm->reverse_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;i<nlocal+nghost;i++)
for (j=0;j<3;j++) {
f[i][j] = 0.0;
torque[i][j] = 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;
// Find angular velocity
wi[0] = omega[i][0];
wi[1] = omega[i][1];
wi[2] = omega[i][2];
// Contribution due to the isotropic terms
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];
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);
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 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
// 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]);
// 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];
// 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_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)*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 = pre[0]*(beta[0][1]*beta[0][1]/(beta[1][1]*beta[1][1]*h_sep));
// Relative velocity at the point of closest approach
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;
// 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;
// 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;
// add to total
torque[i][0] -= vxmu2f*tx;
torque[i][1] -= vxmu2f*ty;
torque[i][2] -= vxmu2f*tz;
}
}
}
}
}
/* ----------------------------------------------------------------------
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 PairLubricateUPoly::compute_RE(double **x)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
int *ilist,*jlist,*numneigh,**firstneigh;
int *type = atom->type;
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;
}
c4science · Help