jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@8126 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2012-05-22 15:08:24 +00:00
parent 8752e626b0
commit 597e68a377
10 changed files with 95 additions and 903 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

133
src/FLD/pair_brownian.cpp
View File

@ -29,24 +29,14 @@
#include "neigh_request.h"
#include "domain.h"
#include "update.h"
#include "modify.h"
#include "fix.h"
#include "fix_deform.h"
#include "fix_wall.h"
#include "input.h"
#include "variable.h"
#include "memory.h"
#include "random_mars.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace MathConst;
// same as fix_wall.cpp
enum{EDGE,CONSTANT,VARIABLE};
/* ---------------------------------------------------------------------- */
PairBrownian::PairBrownian(LAMMPS *lmp) : Pair(lmp)
@ -96,48 +86,7 @@ void PairBrownian::compute(int eflag, int vflag)
double xl[3],a_sq,a_sh,a_pu,Fbmag;
double p1[3],p2[3],p3[3];
int overlaps = 0;
// This section of code adjusts R0/RT0/RS0 if necessary due to changes
// in the volume fraction as a result of fix deform or moving walls
double dims[3], wallcoord;
if (flagVF) // Flag for volume fraction corrections
if (flagdeform || flagwall == 2){ // Possible changes in volume fraction
if (flagdeform && !flagwall)
for (j = 0; j < 3; j++)
dims[j] = domain->prd[j];
else if (flagwall == 2 || (flagdeform && flagwall == 1)){
double wallhi[3], walllo[3];
for (int j = 0; j < 3; j++){
wallhi[j] = domain->prd[j];
walllo[j] = 0;
}
for (int m = 0; m < wallfix->nwall; m++){
int dim = wallfix->wallwhich[m] / 2;
int side = wallfix->wallwhich[m] % 2;
if (wallfix->wallstyle[m] == VARIABLE){
wallcoord = input->variable->compute_equal(wallfix->varindex[m]);
}
else wallcoord = wallfix->coord0[m];
if (side == 0) walllo[dim] = wallcoord;
else wallhi[dim] = wallcoord;
}
for (int j = 0; j < 3; j++)
dims[j] = wallhi[j] - walllo[j];
}
double vol_T = dims[0]*dims[1]*dims[2];
double vol_f = vol_P/vol_T;
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);
}
}
// scale factor for Brownian moments
prethermostat = sqrt(24.0*force->boltz*t_target/update->dt);
@ -171,8 +120,6 @@ void PairBrownian::compute(int eflag, int vflag)
}
}
if (!flagHI) continue;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
@ -371,7 +318,7 @@ void PairBrownian::allocate()
void PairBrownian::settings(int narg, char **arg)
{
if (narg < 7 || narg > 9) error->all(FLERR,"Illegal pair_style command");
if (narg != 7) error->all(FLERR,"Illegal pair_style command");
mu = atof(arg[0]);
flaglog = atoi(arg[1]);
@ -381,17 +328,6 @@ void PairBrownian::settings(int narg, char **arg)
t_target = atof(arg[5]);
seed = atoi(arg[6]);
flagHI = flagVF = 1;
if (narg >= 8) flagHI = atoi(arg[7]);
if (narg == 9) flagVF = atoi(arg[8]);
if (flaglog == 1 && flagHI == 0) {
error->warning(FLERR,"Cannot include log terms without 1/r terms; "
"setting flagHI to 1");
flagHI = 1;
}
// initialize Marsaglia RNG with processor-unique seed
delete random;
@ -476,7 +412,7 @@ void PairBrownian::init_style()
// require monodisperse system with same radii for all types
double radtype;
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");
@ -487,66 +423,17 @@ void PairBrownian::init_style()
// set the isotropic constants that depend on the volume fraction
// vol_T = total volume
// check for fix deform, if exists it must use "remap v"
// If box will change volume, set appropriate flag so that volume
// and v.f. corrections are re-calculated at every step.
//
// If available volume is different from box volume
// due to walls, set volume appropriately; if walls will
// move, set appropriate flag so that volume and v.f. corrections
// are re-calculated at every step.
flagdeform = flagwall = 0;
for (int i = 0; i < modify->nfix; i++){
if (strcmp(modify->fix[i]->style,"deform") == 0)
flagdeform = 1;
else if (strstr(modify->fix[i]->style,"wall") != NULL){
flagwall = 1; // Walls exist
if (((FixWall *) modify->fix[i])->varflag ) {
flagwall = 2; // Moving walls exist
wallfix = (FixWall *) modify->fix[i];
}
}
}
double vol_T = domain->xprd*domain->yprd*domain->zprd;
// set the isotropic constants depending on the volume fraction
// vol_T = total volumeshearing = flagdeform = flagwall = 0;
double vol_T, wallcoord;
if (!flagwall) vol_T = domain->xprd*domain->yprd*domain->zprd;
else {
double wallhi[3], walllo[3];
for (int j = 0; j < 3; j++){
wallhi[j] = domain->prd[j];
walllo[j] = 0;
}
for (int m = 0; m < wallfix->nwall; m++){
int dim = wallfix->wallwhich[m] / 2;
int side = wallfix->wallwhich[m] % 2;
if (wallfix->wallstyle[m] == VARIABLE){
wallfix->varindex[m] = input->variable->find(wallfix->varstr[m]);
// Since fix->wall->init happens after pair->init_style
wallcoord = input->variable->compute_equal(wallfix->varindex[m]);
}
else wallcoord = wallfix->coord0[m];
if (side == 0) walllo[dim] = wallcoord;
else wallhi[dim] = wallcoord;
}
vol_T = (wallhi[0] - walllo[0]) * (wallhi[1] - walllo[1]) *
(wallhi[2] - walllo[2]);
}
// vol_P = volume of particles, assuming mono-dispersity
// vol_f = volume fraction
vol_P = atom->natoms*(4.0/3.0)*MY_PI*pow(rad,3);
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 (!flagVF) vol_f = 0;
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
@ -632,8 +519,6 @@ void PairBrownian::write_restart_settings(FILE *fp)
fwrite(&seed,sizeof(int),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
fwrite(&flagHI,sizeof(int),1,fp);
fwrite(&flagVF,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
@ -653,8 +538,6 @@ void PairBrownian::read_restart_settings(FILE *fp)
fread(&seed, sizeof(int),1,fp);
fread(&offset_flag,sizeof(int),1,fp);
fread(&mix_flag,sizeof(int),1,fp);
fread(&flagHI,sizeof(int),1,fp);
fread(&flagVF,sizeof(int),1,fp);
}
MPI_Bcast(&mu,1,MPI_DOUBLE,0,world);
MPI_Bcast(&flaglog,1,MPI_INT,0,world);
@ -665,8 +548,6 @@ void PairBrownian::read_restart_settings(FILE *fp)
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);
MPI_Bcast(&flagHI,1,MPI_INT,0,world);
MPI_Bcast(&flagVF,1,MPI_INT,0,world);
// additional setup based on restart parameters

6
src/FLD/pair_brownian.h
View File

@ -42,12 +42,6 @@ class PairBrownian : public Pair {
double cut_inner_global,cut_global;
double t_target,mu;
int flaglog,flagfld;
int flagHI, flagVF;
int flagdeform, flagwall;
double vol_P;
double rad;
class FixWall *wallfix;
int seed;
double **cut_inner,**cut;
double R0,RT0;

90
src/FLD/pair_brownian_poly.cpp
View File

@ -30,24 +30,14 @@
#include "neigh_request.h"
#include "domain.h"
#include "update.h"
#include "modify.h"
#include "fix.h"
#include "fix_deform.h"
#include "fix_wall.h"
#include "input.h"
#include "variable.h"
#include "memory.h"
#include "random_mars.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace MathConst;
// same as fix_wall.cpp
enum{EDGE,CONSTANT,VARIABLE};
/* ---------------------------------------------------------------------- */
PairBrownianPoly::PairBrownianPoly(LAMMPS *lmp) : PairBrownian(lmp)
@ -82,47 +72,6 @@ void PairBrownianPoly::compute(int eflag, int vflag)
double xl[3],a_sq,a_sh,a_pu,Fbmag;
double p1[3],p2[3],p3[3];
// This section of code adjusts R0/RT0/RS0 if necessary due to changes
// in the volume fraction as a result of fix deform or moving walls
double dims[3], wallcoord;
if (flagVF) // Flag for volume fraction corrections
if (flagdeform || flagwall == 2){ // Possible changes in volume fraction
if (flagdeform && !flagwall)
for (j = 0; j < 3; j++)
dims[j] = domain->prd[j];
else if (flagwall == 2 || (flagdeform && flagwall == 1)){
double wallhi[3], walllo[3];
for (j = 0; j < 3; j++){
wallhi[j] = domain->prd[j];
walllo[j] = 0;
}
for (int m = 0; m < wallfix->nwall; m++){
int dim = wallfix->wallwhich[m] / 2;
int side = wallfix->wallwhich[m] % 2;
if (wallfix->wallstyle[m] == VARIABLE){
wallcoord = input->variable->compute_equal(wallfix->varindex[m]);
}
else wallcoord = wallfix->coord0[m];
if (side == 0) walllo[dim] = wallcoord;
else wallhi[dim] = wallcoord;
}
for (j = 0; j < 3; j++)
dims[j] = wallhi[j] - walllo[j];
}
double vol_T = dims[0]*dims[1]*dims[2];
double vol_f = vol_P/vol_T;
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);
}
}
// scale factor for Brownian moments
prethermostat = sqrt(24.0*force->boltz*t_target/update->dt);
@ -133,7 +82,6 @@ void PairBrownianPoly::compute(int eflag, int vflag)
numneigh = list->numneigh;
firstneigh = list->firstneigh;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
@ -160,8 +108,6 @@ void PairBrownianPoly::compute(int eflag, int vflag)
}
}
if (!flagHI) continue;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
@ -351,44 +297,20 @@ void PairBrownianPoly::init_style()
// set the isotropic constants that depend on the volume fraction
// vol_T = total volume
double vol_T, wallcoord;
if (!flagwall) vol_T = domain->xprd*domain->yprd*domain->zprd;
else {
double wallhi[3], walllo[3];
for (int j = 0; j < 3; j++){
wallhi[j] = domain->prd[j];
walllo[j] = 0;
}
for (int m = 0; m < wallfix->nwall; m++){
int dim = wallfix->wallwhich[m] / 2;
int side = wallfix->wallwhich[m] % 2;
if (wallfix->wallstyle[m] == VARIABLE){
wallfix->varindex[m] = input->variable->find(wallfix->varstr[m]);
// Since fix->wall->init happens after pair->init_style
wallcoord = input->variable->compute_equal(wallfix->varindex[m]);
}
else wallcoord = wallfix->coord0[m];
if (side == 0) walllo[dim] = wallcoord;
else wallhi[dim] = wallcoord;
}
vol_T = (wallhi[0] - walllo[0]) * (wallhi[1] - walllo[1]) *
(wallhi[2] - walllo[2]);
}
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);
MPI_Allreduce(&volP,&vol_P,1,MPI_DOUBLE,MPI_SUM,world);
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;
if (!flagVF) vol_f = 0;
// set isotropic constants
if (flaglog == 0) {

147
src/FLD/pair_lubricate.cpp
View File

@ -33,12 +33,9 @@
#include "modify.h"
#include "fix.h"
#include "fix_deform.h"
#include "fix_wall.h"
#include "input.h"
#include "variable.h"
#include "memory.h"
#include "random_mars.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
@ -48,10 +45,6 @@ using namespace MathConst;
enum{NO_REMAP,X_REMAP,V_REMAP};
// same as fix_wall.cpp
enum{EDGE,CONSTANT,VARIABLE};
/* ---------------------------------------------------------------------- */
PairLubricate::PairLubricate(LAMMPS *lmp) : Pair(lmp)
@ -164,53 +157,6 @@ void PairLubricate::compute(int eflag, int vflag)
comm->forward_comm_pair(this);
}
// This section of code adjusts R0/RT0/RS0 if necessary due to changes
// in the volume fraction as a result of fix deform or moving walls
double dims[3], wallcoord;
if (flagVF) // Flag for volume fraction corrections
if (flagdeform || flagwall == 2){ // Possible changes in volume fraction
if (flagdeform && !flagwall)
for (j = 0; j < 3; j++)
dims[j] = domain->prd[j];
else if (flagwall == 2 || (flagdeform && flagwall == 1)){
double wallhi[3], walllo[3];
for (int j = 0; j < 3; j++){
wallhi[j] = domain->prd[j];
walllo[j] = 0;
}
for (int m = 0; m < wallfix->nwall; m++){
int dim = wallfix->wallwhich[m] / 2;
int side = wallfix->wallwhich[m] % 2;
if (wallfix->wallstyle[m] == VARIABLE){
wallcoord = input->variable->compute_equal(wallfix->varindex[m]);
}
else wallcoord = wallfix->coord0[m];
if (side == 0) walllo[dim] = wallcoord;
else wallhi[dim] = wallcoord;
}
for (int j = 0; j < 3; j++)
dims[j] = wallhi[j] - walllo[j];
}
double vol_T = dims[0]*dims[1]*dims[2];
double vol_f = vol_P/vol_T;
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);
}
}
// end of R0 adjustment code
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
@ -246,8 +192,6 @@ void PairLubricate::compute(int eflag, int vflag)
}
}
if (!flagHI) continue;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
@ -487,7 +431,7 @@ void PairLubricate::allocate()
void PairLubricate::settings(int narg, char **arg)
{
if (narg < 5 || narg > 7) error->all(FLERR,"Illegal pair_style command");
if (narg != 5) error->all(FLERR,"Illegal pair_style command");
mu = atof(arg[0]);
flaglog = atoi(arg[1]);
@ -495,18 +439,6 @@ void PairLubricate::settings(int narg, char **arg)
cut_inner_global = atof(arg[3]);
cut_global = atof(arg[4]);
flagVF = flagHI = 1;
if (narg >= 6) flagHI = atoi(arg[5]);
if (narg == 7) flagVF = atoi(arg[6]);
if (flaglog == 1 && flagHI == 0) {
error->warning(FLERR,"Cannot include log terms without 1/r terms; "
"setting flagHI to 1");
flagHI = 1;
}
// reset cutoffs that have been explicitly set
if (allocated) {
@ -570,7 +502,7 @@ void PairLubricate::init_style()
// require that atom radii are identical within each type
// require monodisperse system with same radii for all types
double radtype;
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");
@ -579,69 +511,16 @@ void PairLubricate::init_style()
rad = radtype;
}
// check for fix deform, if exists it must use "remap v"
// If box will change volume, set appropriate flag so that volume
// and v.f. corrections are re-calculated at every step.
//
// If available volume is different from box volume
// due to walls, set volume appropriately; if walls will
// move, set appropriate flag so that volume and v.f. corrections
// are re-calculated at every step.
shearing = flagdeform = flagwall = 0;
for (int i = 0; i < modify->nfix; i++){
if (strcmp(modify->fix[i]->style,"deform") == 0) {
shearing = flagdeform = 1;
if (((FixDeform *) modify->fix[i])->remapflag != V_REMAP)
error->all(FLERR,"Using pair lubricate with inconsistent "
"fix deform remap option");
}
if (strstr(modify->fix[i]->style,"wall") != NULL){
flagwall = 1; // Walls exist
if (((FixWall *) modify->fix[i])->varflag ) {
flagwall = 2; // Moving walls exist
wallfix = (FixWall *) modify->fix[i];
}
}
}
// set the isotropic constants that depend on the volume fraction
// vol_T = total volume
double vol_T;
double wallcoord;
if (!flagwall) vol_T = domain->xprd*domain->yprd*domain->zprd;
else {
double wallhi[3], walllo[3];
for (int j = 0; j < 3; j++){
wallhi[j] = domain->prd[j];
walllo[j] = 0;
}
for (int m = 0; m < wallfix->nwall; m++){
int dim = wallfix->wallwhich[m] / 2;
int side = wallfix->wallwhich[m] % 2;
if (wallfix->wallstyle[m] == VARIABLE){
wallfix->varindex[m] = input->variable->find(wallfix->varstr[m]);
//Since fix->wall->init happens after pair->init_style
wallcoord = input->variable->compute_equal(wallfix->varindex[m]);
}
else wallcoord = wallfix->coord0[m];
if (side == 0) walllo[dim] = wallcoord;
else wallhi[dim] = wallcoord;
}
vol_T = (wallhi[0] - walllo[0]) * (wallhi[1] - walllo[1]) *
(wallhi[2] - walllo[2]);
}
double vol_T = domain->xprd*domain->yprd*domain->zprd;
// vol_P = volume of particles, assuming monodispersity
// vol_f = volume fraction
vol_P = atom->natoms*(4.0/3.0)*MY_PI*pow(rad,3);
double vol_P = atom->natoms*(4.0/3.0)*MY_PI*pow(rad,3);
double vol_f = vol_P/vol_T;
if (!flagVF) vol_f = 0;
// set isotropic constants for FLD
@ -655,6 +534,16 @@ void PairLubricate::init_style()
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
@ -737,8 +626,6 @@ void PairLubricate::write_restart_settings(FILE *fp)
fwrite(&cut_global,sizeof(double),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
fwrite(&flagHI,sizeof(int),1,fp);
fwrite(&flagVF,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
@ -756,8 +643,6 @@ void PairLubricate::read_restart_settings(FILE *fp)
fread(&cut_global,sizeof(double),1,fp);
fread(&offset_flag,sizeof(int),1,fp);
fread(&mix_flag,sizeof(int),1,fp);
fread(&flagHI,sizeof(int),1,fp);
fread(&flagVF,sizeof(int),1,fp);
}
MPI_Bcast(&mu,1,MPI_DOUBLE,0,world);
MPI_Bcast(&flaglog,1,MPI_INT,0,world);
@ -766,8 +651,6 @@ void PairLubricate::read_restart_settings(FILE *fp)
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);
MPI_Bcast(&flagHI,1,MPI_INT,0,world);
MPI_Bcast(&flagVF,1,MPI_INT,0,world);
}
/* ---------------------------------------------------------------------- */

6
src/FLD/pair_lubricate.h
View File

@ -45,13 +45,7 @@ class PairLubricate : public Pair {
protected:
double mu,cut_inner_global,cut_global;
double rad;
int flaglog,flagfld,shearing;
int flagdeform, flagwall;
double vol_P;
class FixWall *wallfix;
int flagVF, flagHI;
double Ef[3][3];
double R0,RT0,RS0;
double **cut_inner,**cut;

245
src/FLD/pair_lubricateU.cpp
View File

@ -31,25 +31,14 @@
#include "domain.h"
#include "update.h"
#include "math_const.h"
#include "modify.h"
#include "fix.h"
#include "fix_deform.h"
#include "fix_wall.h"
#include "input.h"
#include "variable.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace MathConst;
#define TOL 1E-4 // tolerance for conjugate gradient
// same as fix_wall.cpp
enum{EDGE,CONSTANT,VARIABLE};
/* ---------------------------------------------------------------------- */
PairLubricateU::PairLubricateU(LAMMPS *lmp) : Pair(lmp)
@ -609,53 +598,8 @@ void PairLubricateU::compute_Fh(double **x)
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// This section of code adjusts R0/RT0/RS0 if necessary due to changes
// in the volume fraction as a result of fix deform or moving walls
double dims[3], wallcoord;
if (flagVF) // Flag for volume fraction corrections
if (flagdeform || flagwall == 2){ // Possible changes in volume fraction
if (flagdeform && !flagwall)
for (j = 0; j < 3; j++)
dims[j] = domain->prd[j];
else if (flagwall == 2 || (flagdeform && flagwall == 1)){
double wallhi[3], walllo[3];
for (int j = 0; j < 3; j++){
wallhi[j] = domain->prd[j];
walllo[j] = 0;
}
for (int m = 0; m < wallfix->nwall; m++){
int dim = wallfix->wallwhich[m] / 2;
int side = wallfix->wallwhich[m] % 2;
if (wallfix->wallstyle[m] == VARIABLE){
wallcoord = input->variable->compute_equal(wallfix->varindex[m]);
}
else wallcoord = wallfix->coord0[m];
if (side == 0) walllo[dim] = wallcoord;
else wallhi[dim] = wallcoord;
}
for (int j = 0; j < 3; j++)
dims[j] = wallhi[j] - walllo[j];
}
double vol_T = dims[0]*dims[1]*dims[2];
double vol_f = vol_P/vol_T;
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);
}
}
// end of R0 adjustment code
// 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;
@ -694,8 +638,6 @@ void PairLubricateU::compute_Fh(double **x)
wi[1] = omega[i][1];
wi[2] = omega[i][2];
if (!flagHI) continue;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
@ -848,50 +790,8 @@ void PairLubricateU::compute_RU()
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// This section of code adjusts R0/RT0/RS0 if necessary due to changes
// in the volume fraction as a result of fix deform or moving walls
double dims[3], wallcoord;
if (flagVF) // Flag for volume fraction corrections
if (flagdeform || flagwall == 2){ // Possible changes in volume fraction
if (flagdeform && !flagwall)
for (j = 0; j < 3; j++)
dims[j] = domain->prd[j];
else if (flagwall == 2 || (flagdeform && flagwall == 1)){
double wallhi[3], walllo[3];
for (int j = 0; j < 3; j++){
wallhi[j] = domain->prd[j];
walllo[j] = 0;
}
for (int m = 0; m < wallfix->nwall; m++){
int dim = wallfix->wallwhich[m] / 2;
int side = wallfix->wallwhich[m] % 2;
if (wallfix->wallstyle[m] == VARIABLE){
wallcoord = input->variable->compute_equal(wallfix->varindex[m]);
}
else wallcoord = wallfix->coord0[m];
if (side == 0) walllo[dim] = wallcoord;
else wallhi[dim] = wallcoord;
}
for (int j = 0; j < 3; j++)
dims[j] = wallhi[j] - walllo[j];
}
double vol_T = dims[0]*dims[1]*dims[2];
double vol_f = vol_P/vol_T;
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);
}
}
// end of R0 adjustment code
// Initialize f to zero
for (i=0;i<nlocal+nghost;i++)
for (j=0;j<3;j++) {
f[i][j] = 0.0;
@ -924,8 +824,6 @@ void PairLubricateU::compute_RU()
torque[i][1] += -vxmu2f*RT0*wi[1];
torque[i][2] += -vxmu2f*RT0*wi[2];
if (!flagHI) continue;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
@ -1125,50 +1023,8 @@ void PairLubricateU::compute_RU(double **x)
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// This section of code adjusts R0/RT0/RS0 if necessary due to changes
// in the volume fraction as a result of fix deform or moving walls
double dims[3], wallcoord;
if (flagVF) // Flag for volume fraction corrections
if (flagdeform || flagwall == 2){ // Possible changes in volume fraction
if (flagdeform && !flagwall)
for (j = 0; j < 3; j++)
dims[j] = domain->prd[j];
else if (flagwall == 2 || (flagdeform && flagwall == 1)){
double wallhi[3], walllo[3];
for (int j = 0; j < 3; j++){
wallhi[j] = domain->prd[j];
walllo[j] = 0;
}
for (int m = 0; m < wallfix->nwall; m++){
int dim = wallfix->wallwhich[m] / 2;
int side = wallfix->wallwhich[m] % 2;
if (wallfix->wallstyle[m] == VARIABLE){
wallcoord = input->variable->compute_equal(wallfix->varindex[m]);
}
else wallcoord = wallfix->coord0[m];
if (side == 0) walllo[dim] = wallcoord;
else wallhi[dim] = wallcoord;
}
for (int j = 0; j < 3; j++)
dims[j] = wallhi[j] - walllo[j];
}
double vol_T = dims[0]*dims[1]*dims[2];
double vol_f = vol_P/vol_T;
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);
}
}
// end of R0 adjustment code
// Initialize f to zero
for (i=0;i<nlocal+nghost;i++)
for (j=0;j<3;j++) {
f[i][j] = 0.0;
@ -1201,8 +1057,6 @@ void PairLubricateU::compute_RU(double **x)
torque[i][1] += -vxmu2f*RT0*wi[1];
torque[i][2] += -vxmu2f*RT0*wi[2];
if (!flagHI) continue;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
@ -1406,8 +1260,6 @@ void PairLubricateU::compute_RE()
numneigh = list->numneigh;
firstneigh = list->firstneigh;
if (!flagHI) return;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
@ -1582,8 +1434,6 @@ void PairLubricateU::compute_RE(double **x)
int overlaps = 0;
double vi[3],vj[3],wi[3],wj[3],xl[3],a_sq,a_sh,a_pu,Fbmag,del,delmin,eta;
if (!flagHI) return;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
@ -1759,22 +1609,14 @@ void PairLubricateU::allocate()
void PairLubricateU::settings(int narg, char **arg)
{
if (narg < 5 || narg > 7) error->all(FLERR,"Illegal pair_style command");
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]);
flagHI = flagVF = 1;
if (narg >= 6) flagHI = atoi(arg[5]);
if (narg == 7) flagVF = atoi(arg[6]);
if (flaglog == 1 && flagHI == 0) {
error->warning(FLERR,"Cannot include log terms without 1/r terms. Setting flagHI to 1.");
flagHI = 1;
}
// reset cutoffs that have been explicitly set
if (allocated) {
@ -1851,90 +1693,43 @@ void PairLubricateU::init_style()
// require that atom radii are identical within each type
// require monodisperse system with same radii for all types
double radi, radtype;
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");
radi = radtype;
}
// check for fix deform, if exists it must use "remap v"
// If box will change volume, set appropriate flag so that volume
// and v.f. corrections are re-calculated at every step.
//
// If available volume is different from box volume
// due to walls, set volume appropriately; if walls will
// move, set appropriate flag so that volume and v.f. corrections
// are re-calculated at every step.
flagdeform = flagwall = 0;
for (int i = 0; i < modify->nfix; i++){
if (strcmp(modify->fix[i]->style,"deform") == 0)
flagdeform = 1;
else if (strstr(modify->fix[i]->style,"wall") != NULL){
flagwall = 1; // Walls exist
if (((FixWall *) modify->fix[i])->varflag ) {
flagwall = 2; // Moving walls exist
wallfix = (FixWall *) modify->fix[i];
}
}
rad = radtype;
}
// set the isotropic constants depending on the volume fraction
// vol_T = total volumeshearing = flagdeform = flagwall = 0;
double vol_T, wallcoord;
if (!flagwall) vol_T = domain->xprd*domain->yprd*domain->zprd;
else {
double wallhi[3], walllo[3];
for (int j = 0; j < 3; j++){
wallhi[j] = domain->prd[j];
walllo[j] = 0;
}
for (int m = 0; m < wallfix->nwall; m++){
int dim = wallfix->wallwhich[m] / 2;
int side = wallfix->wallwhich[m] % 2;
if (wallfix->wallstyle[m] == VARIABLE){
wallfix->varindex[m] = input->variable->find(wallfix->varstr[m]);
//Since fix->wall->init happens after pair->init_style
wallcoord = input->variable->compute_equal(wallfix->varindex[m]);
}
// vol_T = total volume
else wallcoord = wallfix->coord0[m];
if (side == 0) walllo[dim] = wallcoord;
else wallhi[dim] = wallcoord;
}
vol_T = (wallhi[0] - walllo[0]) * (wallhi[1] - walllo[1]) *
(wallhi[2] - walllo[2]);
}
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];
MPI_Allreduce(&tmp,&rad,1,MPI_DOUBLE,MPI_MAX,world);
double radi;
MPI_Allreduce(&tmp,&radi,1,MPI_DOUBLE,MPI_MAX,world);
vol_P = atom->natoms * (4.0/3.0)*MY_PI*pow(rad,3);
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;
if (!flagVF) vol_f = 0;
// set the isotropic constant
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
RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f);
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*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);
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);
}
}
@ -2011,8 +1806,6 @@ void PairLubricateU::write_restart_settings(FILE *fp)
fwrite(&cut_global,sizeof(double),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
fwrite(&flagHI,sizeof(int),1,fp);
fwrite(&flagVF,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
@ -2029,8 +1822,6 @@ void PairLubricateU::read_restart_settings(FILE *fp)
fread(&cut_global,sizeof(double),1,fp);
fread(&offset_flag,sizeof(int),1,fp);
fread(&mix_flag,sizeof(int),1,fp);
fread(&flagHI,sizeof(int),1,fp);
fread(&flagVF,sizeof(int),1,fp);
}
MPI_Bcast(&mu,1,MPI_DOUBLE,0,world);
MPI_Bcast(&flaglog,1,MPI_INT,0,world);
@ -2038,8 +1829,6 @@ void PairLubricateU::read_restart_settings(FILE *fp)
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);
MPI_Bcast(&flagHI,1,MPI_INT,0,world);
MPI_Bcast(&flagVF,1,MPI_INT,0,world);
}
/*---------------------------------------------------------------------------*/

6
src/FLD/pair_lubricateU.h
View File

@ -43,13 +43,7 @@ class PairLubricateU : public Pair {
protected:
double cut_inner_global,cut_global;
double mu;
double rad;
int flaglog;
int flagdeform, flagwall;
int flagVF, flagHI;
double vol_P;
class FixWall *wallfix;
double gdot,Ef[3][3];
double **cut_inner,**cut;
void allocate();

249
src/FLD/pair_lubricateU_poly.cpp
View File

@ -32,12 +32,6 @@
#include "neigh_request.h"
#include "domain.h"
#include "update.h"
#include "modify.h"
#include "fix.h"
#include "fix_deform.h"
#include "fix_wall.h"
#include "input.h"
#include "variable.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
@ -47,11 +41,6 @@ using namespace MathConst;
#define TOL 1E-3 // tolerance for conjugate gradient
// same as fix_wall.cpp
enum{EDGE,CONSTANT,VARIABLE};
/* ---------------------------------------------------------------------- */
PairLubricateUPoly::PairLubricateUPoly(LAMMPS *lmp) :
@ -361,48 +350,6 @@ void PairLubricateUPoly::compute_Fh(double **x)
beta[0][0] = beta[1][0] = beta[1][4] = 0.0;
// This section of code adjusts R0/RT0/RS0 if necessary due to changes
// in the volume fraction as a result of fix deform or moving walls
double dims[3], wallcoord;
if (flagVF) // Flag for volume fraction corrections
if (flagdeform || flagwall == 2){ // Possible changes in volume fraction
if (flagdeform && !flagwall)
for (j = 0; j < 3; j++)
dims[j] = domain->prd[j];
else if (flagwall == 2 || (flagdeform && flagwall == 1)){
double wallhi[3], walllo[3];
for (int j = 0; j < 3; j++){
wallhi[j] = domain->prd[j];
walllo[j] = 0;
}
for (int m = 0; m < wallfix->nwall; m++){
int dim = wallfix->wallwhich[m] / 2;
int side = wallfix->wallwhich[m] % 2;
if (wallfix->wallstyle[m] == VARIABLE){
wallcoord = input->variable->compute_equal(wallfix->varindex[m]);
}
else wallcoord = wallfix->coord0[m];
if (side == 0) walllo[dim] = wallcoord;
else wallhi[dim] = wallcoord;
}
for (int j = 0; j < 3; j++)
dims[j] = wallhi[j] - walllo[j];
}
double vol_T = dims[0]*dims[1]*dims[2];
double vol_f = vol_P/vol_T;
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);
}
}
// end of R0 adjustment code
// Set force to zero which is the final value after this pair interaction
for (i=0;i<nlocal+nghost;i++)
@ -426,7 +373,7 @@ void PairLubricateUPoly::compute_Fh(double **x)
radi = radius[i];
jlist = firstneigh[i];
jnum = numneigh[i];
pre[1] = 8.0*(pre[0] = MY_PI*mu*radi)*radi*radi; // BROKEN?? Should be "+"??
pre[1] = 8.0*(pre[0] = MY_PI*mu*radi)*radi*radi;
pre[0] *= 6.0;
// Find the contribution to stress from isotropic RS0
@ -642,49 +589,6 @@ void PairLubricateUPoly::compute_RU(double **x)
beta[0][0] = beta[1][0] = beta[1][4] = 0.0;
// This section of code adjusts R0/RT0/RS0 if necessary due to changes
// in the volume fraction as a result of fix deform or moving walls
double dims[3], wallcoord;
if (flagVF) // Flag for volume fraction corrections
if (flagdeform || flagwall == 2){ // Possible changes in volume fraction
if (flagdeform && !flagwall)
for (j = 0; j < 3; j++)
dims[j] = domain->prd[j];
else if (flagwall == 2 || (flagdeform && flagwall == 1)){
double wallhi[3], walllo[3];
for (j = 0; j < 3; j++){
wallhi[j] = domain->prd[j];
walllo[j] = 0;
}
for (int m = 0; m < wallfix->nwall; m++){
int dim = wallfix->wallwhich[m] / 2;
int side = wallfix->wallwhich[m] % 2;
if (wallfix->wallstyle[m] == VARIABLE){
wallcoord = input->variable->compute_equal(wallfix->varindex[m]);
}
else wallcoord = wallfix->coord0[m];
if (side == 0) walllo[dim] = wallcoord;
else wallhi[dim] = wallcoord;
}
for (j = 0; j < 3; j++)
dims[j] = wallhi[j] - walllo[j];
}
double vol_T = dims[0]*dims[1]*dims[2];
double vol_f = vol_P/vol_T;
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);
}
}
// end of R0 adjustment code
// Initialize f to zero
for (i=0;i<nlocal+nghost;i++)
@ -720,8 +624,7 @@ void PairLubricateUPoly::compute_RU(double **x)
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 (!flagHI) continue;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
@ -938,9 +841,7 @@ void PairLubricateUPoly::compute_RE(double **x)
double a_sq = 0.0;
double a_sh = 0.0;
double a_pu = 0.0;
if (!flagHI) return;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
@ -961,6 +862,7 @@ void PairLubricateUPoly::compute_RE(double **x)
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];
@ -1122,9 +1024,10 @@ void PairLubricateUPoly::compute_RE(double **x)
void PairLubricateUPoly::settings(int narg, char **arg)
{
double vol_T;
int itype;
if (narg < 5 || narg > 7) error->all(FLERR,"Illegal pair_style command");
if (narg != 5) error->all(FLERR,"Illegal pair_style command");
mu = atof(arg[0]);
flaglog = atoi(arg[1]);
@ -1132,16 +1035,6 @@ void PairLubricateUPoly::settings(int narg, char **arg)
cut_global = atof(arg[3]);
gdot = atof(arg[4]);
flagHI = flagVF = 1;
if (narg >= 6) flagHI = atoi(arg[5]);
if (narg == 7) flagVF = atoi(arg[6]);
if (flaglog == 1 && flagHI == 0) {
error->warning(FLERR,"Cannot include log terms without 1/r terms; "
"setting flagHI to 1");
flagHI = 1;
}
// reset cutoffs that have been explicitly set
if (allocated) {
@ -1166,7 +1059,48 @@ void PairLubricateUPoly::settings(int narg, char **arg)
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);
}
}
/* ----------------------------------------------------------------------
@ -1193,95 +1127,6 @@ void PairLubricateUPoly::init_style()
for (int i = 0; i < nlocal; i++)
if (radius[i] == 0.0)
error->one(FLERR,"Pair lubricate/poly requires extended particles");
// Set the isotropic constants depending on the volume fraction
// Find the total volume
// check for fix deform, if exists it must use "remap v"
// If box will change volume, set appropriate flag so that volume
// and v.f. corrections are re-calculated at every step.
//
// If available volume is different from box volume
// due to walls, set volume appropriately; if walls will
// move, set appropriate flag so that volume and v.f. corrections
// are re-calculated at every step.
flagdeform = flagwall = 0;
for (int i = 0; i < modify->nfix; i++){
if (strcmp(modify->fix[i]->style,"deform") == 0)
flagdeform = 1;
else if (strstr(modify->fix[i]->style,"wall") != NULL){
flagwall = 1; // Walls exist
if (((FixWall *) modify->fix[i])->varflag ) {
flagwall = 2; // Moving walls exist
wallfix = (FixWall *) modify->fix[i];
}
}
}
// set the isotropic constants depending on the volume fraction
// vol_T = total volumeshearing = flagdeform = flagwall = 0;
double vol_T, wallcoord;
if (!flagwall) vol_T = domain->xprd*domain->yprd*domain->zprd;
else {
double wallhi[3], walllo[3];
for (int j = 0; j < 3; j++){
wallhi[j] = domain->prd[j];
walllo[j] = 0;
}
for (int m = 0; m < wallfix->nwall; m++){
int dim = wallfix->wallwhich[m] / 2;
int side = wallfix->wallwhich[m] % 2;
if (wallfix->wallstyle[m] == VARIABLE){
wallfix->varindex[m] = input->variable->find(wallfix->varstr[m]);
//Since fix->wall->init happens after pair->init_style
wallcoord = input->variable->compute_equal(wallfix->varindex[m]);
}
else wallcoord = wallfix->coord0[m];
if (side == 0) walllo[dim] = wallcoord;
else wallhi[dim] = wallcoord;
}
vol_T = (wallhi[0] - walllo[0]) * (wallhi[1] - walllo[1]) *
(wallhi[2] - walllo[2]);
}
// Assuming monodisperse spheres, find the volume of the particles
double volP = 0.0;
for (int i = 0; i < nlocal; i++)
volP += (4.0/3.0)*MY_PI*pow(atom->radius[i],3);
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 (!flagVF) vol_f = 0;
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);
}
int irequest = neighbor->request(this);
neighbor->requests[irequest]->half = 0;

4
src/FLD/pair_lubricateU_poly.h
View File

@ -33,10 +33,6 @@ class PairLubricateUPoly : public PairLubricateU {
void init_style();
private:
double vol_P;
int flagdeform, flagwall, flagVF, flagHI;
class FixWall *wallfix;
void iterate(double **, int);
void compute_RE(double **);
void compute_RU(double **);

112
src/FLD/pair_lubricate_poly.cpp
View File

@ -36,9 +36,6 @@
#include "fix_deform.h"
#include "memory.h"
#include "random_mars.h"
#include "fix_wall.h"
#include "input.h"
#include "variable.h"
#include "math_const.h"
#include "error.h"
@ -49,11 +46,6 @@ using namespace MathConst;
enum{NO_REMAP,X_REMAP,V_REMAP};
// same as fix_wall.cpp
enum{EDGE,CONSTANT,VARIABLE};
/* ---------------------------------------------------------------------- */
PairLubricatePoly::PairLubricatePoly(LAMMPS *lmp) : PairLubricate(lmp)
@ -148,50 +140,6 @@ void PairLubricatePoly::compute(int eflag, int vflag)
comm->forward_comm_pair(this);
}
// This section of code adjusts R0/RT0/RS0 if necessary due to changes
// in the volume fraction as a result of fix deform or moving walls
double dims[3], wallcoord;
if (flagVF) // Flag for volume fraction corrections
if (flagdeform || flagwall == 2){ // Possible changes in volume fraction
if (flagdeform && !flagwall)
for (j = 0; j < 3; j++)
dims[j] = domain->prd[j];
else if (flagwall == 2 || (flagdeform && flagwall == 1)){
double wallhi[3], walllo[3];
for (int j = 0; j < 3; j++){
wallhi[j] = domain->prd[j];
walllo[j] = 0;
}
for (int m = 0; m < wallfix->nwall; m++){
int dim = wallfix->wallwhich[m] / 2;
int side = wallfix->wallwhich[m] % 2;
if (wallfix->wallstyle[m] == VARIABLE){
wallcoord = input->variable->compute_equal(wallfix->varindex[m]);
}
else wallcoord = wallfix->coord0[m];
if (side == 0) walllo[dim] = wallcoord;
else wallhi[dim] = wallcoord;
}
for (int j = 0; j < 3; j++)
dims[j] = wallhi[j] - walllo[j];
}
double vol_T = dims[0]*dims[1]*dims[2];
double vol_f = vol_P/vol_T;
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);
}
}
// end of R0 adjustment code
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
@ -227,8 +175,6 @@ void PairLubricatePoly::compute(int eflag, int vflag)
}
}
if (!flagHI) continue;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
delx = xtmp - x[j][0];
@ -455,7 +401,7 @@ void PairLubricatePoly::init_style()
if (!atom->sphere_flag)
error->all(FLERR,"Pair lubricate/poly requires atom style sphere");
// ensure all particles are finite-size
// insure all particles are finite-size
// for pair hybrid, should limit test to types using the pair style
double *radius = atom->radius;
@ -473,66 +419,14 @@ void PairLubricatePoly::init_style()
// set the isotropic constants that depend on the volume fraction
// vol_T = total volume
// check for fix deform, if exists it must use "remap v"
// If box will change volume, set appropriate flag so that volume
// and v.f. corrections are re-calculated at every step.
//
// If available volume is different from box volume
// due to walls, set volume appropriately; if walls will
// move, set appropriate flag so that volume and v.f. corrections
// are re-calculated at every step.
shearing = flagdeform = flagwall = 0;
for (int i = 0; i < modify->nfix; i++){
if (strcmp(modify->fix[i]->style,"deform") == 0) {
shearing = flagdeform = 1;
if (((FixDeform *) modify->fix[i])->remapflag != V_REMAP)
error->all(FLERR,"Using pair lubricate with inconsistent "
"fix deform remap option");
}
if (strstr(modify->fix[i]->style,"wall") != NULL){
flagwall = 1; // Walls exist
if (((FixWall *) modify->fix[i])->varflag ) {
flagwall = 2; // Moving walls exist
wallfix = (FixWall *) modify->fix[i];
}
}
}
double vol_T;
double wallcoord;
if (!flagwall) vol_T = domain->xprd*domain->yprd*domain->zprd;
else {
double wallhi[3], walllo[3];
for (int j = 0; j < 3; j++){
wallhi[j] = domain->prd[j];
walllo[j] = 0;
}
for (int m = 0; m < wallfix->nwall; m++){
int dim = wallfix->wallwhich[m] / 2;
int side = wallfix->wallwhich[m] % 2;
if (wallfix->wallstyle[m] == VARIABLE){
wallfix->varindex[m] = input->variable->find(wallfix->varstr[m]);
//Since fix->wall->init happens after pair->init_style
wallcoord = input->variable->compute_equal(wallfix->varindex[m]);
}
else wallcoord = wallfix->coord0[m];
if (side == 0) walllo[dim] = wallcoord;
else wallhi[dim] = wallcoord;
}
vol_T = (wallhi[0] - walllo[0]) * (wallhi[1] - walllo[1]) *
(wallhi[2] - walllo[2]);
}
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;
if (!flagVF) vol_f = 0;
// set isotropic constants