ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add files via upload

Browse Source
This commit is contained in:
To Quy-Dong
2019-11-08 10:51:24 +01:00
committed by GitHub
parent afc9627506
commit f4ff35f2ce
1 changed files with 250 additions and 240 deletions
Download Patch File Download Diff File Expand all files Collapse all files

490
src/fix_wall_reflect.cpp
View File

@ -1,240 +1,250 @@
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under certain rights in this software. This software is distributed under
the GNU General Public License. the GNU General Public License.
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
#include "fix_wall_reflect.h" #include "fix_wall_reflect.h"
#include <cstring> #include <cstring>
#include "atom.h" #include "atom.h"
#include "comm.h" #include "comm.h"
#include "update.h" #include "update.h"
#include "modify.h" #include "modify.h"
#include "domain.h" #include "domain.h"
#include "lattice.h" #include "lattice.h"
#include "input.h" #include "input.h"
#include "variable.h" #include "variable.h"
#include "error.h" #include "error.h"
#include "force.h" #include "force.h"
using namespace LAMMPS_NS;
using namespace FixConst; using namespace LAMMPS_NS;
using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), /* ---------------------------------------------------------------------- */
nwall(0)
{ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) :
if (narg < 4) error->all(FLERR,"Illegal fix wall/reflect command"); Fix(lmp, narg, arg),
nwall(0)
if (strcmp(arg[2],"wall/reflect") == 0) { // child class can be stochastic {
dynamic_group_allow = 1;
if (narg < 4) error->all(FLERR,"Illegal fix wall/reflect command");
// parse args
if (strcmp(arg[2],"wall/reflect/stochastic") == 0) return;// child class can be stochastic
nwall = 0;
int scaleflag = 1;
dynamic_group_allow = 1;
int iarg = 3;
while (iarg < narg) { // parse args
if ((strcmp(arg[iarg],"xlo") == 0) || (strcmp(arg[iarg],"xhi") == 0) ||
(strcmp(arg[iarg],"ylo") == 0) || (strcmp(arg[iarg],"yhi") == 0) || nwall = 0;
(strcmp(arg[iarg],"zlo") == 0) || (strcmp(arg[iarg],"zhi") == 0)) { int scaleflag = 1;
if (iarg+2 > narg)
error->all(FLERR,"Illegal fix wall/reflect command"); int iarg = 3;
while (iarg < narg) {
int newwall; if ((strcmp(arg[iarg],"xlo") == 0) || (strcmp(arg[iarg],"xhi") == 0) ||
if (strcmp(arg[iarg],"xlo") == 0) newwall = XLO; (strcmp(arg[iarg],"ylo") == 0) || (strcmp(arg[iarg],"yhi") == 0) ||
else if (strcmp(arg[iarg],"xhi") == 0) newwall = XHI; (strcmp(arg[iarg],"zlo") == 0) || (strcmp(arg[iarg],"zhi") == 0)) {
else if (strcmp(arg[iarg],"ylo") == 0) newwall = YLO; if (iarg+2 > narg) error->all(FLERR,"Illegal fix wall/reflect command");
else if (strcmp(arg[iarg],"yhi") == 0) newwall = YHI;
else if (strcmp(arg[iarg],"zlo") == 0) newwall = ZLO; int newwall;
else if (strcmp(arg[iarg],"zhi") == 0) newwall = ZHI; if (strcmp(arg[iarg],"xlo") == 0) newwall = XLO;
else if (strcmp(arg[iarg],"xhi") == 0) newwall = XHI;
for (int m = 0; (m < nwall) && (m < 6); m++) else if (strcmp(arg[iarg],"ylo") == 0) newwall = YLO;
if (newwall == wallwhich[m]) else if (strcmp(arg[iarg],"yhi") == 0) newwall = YHI;
error->all(FLERR,"Wall defined twice in fix wall/reflect command"); else if (strcmp(arg[iarg],"zlo") == 0) newwall = ZLO;
else if (strcmp(arg[iarg],"zhi") == 0) newwall = ZHI;
wallwhich[nwall] = newwall;
if (strcmp(arg[iarg+1],"EDGE") == 0) { for (int m = 0; (m < nwall) && (m < 6); m++)
wallstyle[nwall] = EDGE; if (newwall == wallwhich[m])
int dim = wallwhich[nwall] / 2; error->all(FLERR,"Wall defined twice in fix wall/reflect command");
int side = wallwhich[nwall] % 2;
if (side == 0) coord0[nwall] = domain->boxlo[dim]; wallwhich[nwall] = newwall;
else coord0[nwall] = domain->boxhi[dim]; if (strcmp(arg[iarg+1],"EDGE") == 0) {
} else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) { wallstyle[nwall] = EDGE;
wallstyle[nwall] = VARIABLE; int dim = wallwhich[nwall] / 2;
int n = strlen(&arg[iarg+1][2]) + 1; int side = wallwhich[nwall] % 2;
varstr[nwall] = new char[n]; if (side == 0) coord0[nwall] = domain->boxlo[dim];
strcpy(varstr[nwall],&arg[iarg+1][2]); else coord0[nwall] = domain->boxhi[dim];
} else { } else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
wallstyle[nwall] = CONSTANT; wallstyle[nwall] = VARIABLE;
coord0[nwall] = force->numeric(FLERR,arg[iarg+1]); int n = strlen(&arg[iarg+1][2]) + 1;
} varstr[nwall] = new char[n];
strcpy(varstr[nwall],&arg[iarg+1][2]);
nwall++; } else {
iarg += 2; wallstyle[nwall] = CONSTANT;
coord0[nwall] = force->numeric(FLERR,arg[iarg+1]);
} else if (strcmp(arg[iarg],"units") == 0) { }
if (iarg+2 > narg) error->all(FLERR,"Illegal wall/reflect command");
if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; nwall++;
else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; iarg += 2;
else error->all(FLERR,"Illegal fix wall/reflect command");
iarg += 2;
} else error->all(FLERR,"Illegal fix wall/reflect command"); } else if (strcmp(arg[iarg],"units") == 0) {
} if (iarg+2 > narg) error->all(FLERR,"Illegal wall/reflect command");
if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
// error check else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
else error->all(FLERR,"Illegal fix wall/reflect command");
if (nwall == 0) error->all(FLERR,"Illegal fix wall command"); iarg += 2;
} else error->all(FLERR,"Illegal fix wall/reflect command");
for (int m = 0; m < nwall; m++) { }
if ((wallwhich[m] == XLO || wallwhich[m] == XHI) && domain->xperiodic)
error->all(FLERR,"Cannot use fix wall/reflect in periodic dimension"); // error check
if ((wallwhich[m] == YLO || wallwhich[m] == YHI) && domain->yperiodic)
error->all(FLERR,"Cannot use fix wall/reflect in periodic dimension"); if (nwall == 0) error->all(FLERR,"Illegal fix wall command");
if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->zperiodic)
error->all(FLERR,"Cannot use fix wall/reflect in periodic dimension"); for (int m = 0; m < nwall; m++) {
} if ((wallwhich[m] == XLO || wallwhich[m] == XHI) && domain->xperiodic)
error->all(FLERR,"Cannot use fix wall/reflect in periodic dimension");
for (int m = 0; m < nwall; m++) if ((wallwhich[m] == YLO || wallwhich[m] == YHI) && domain->yperiodic)
if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->dimension == 2) error->all(FLERR,"Cannot use fix wall/reflect in periodic dimension");
error->all(FLERR, if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->zperiodic)
"Cannot use fix wall/reflect zlo/zhi for a 2d simulation"); error->all(FLERR,"Cannot use fix wall/reflect in periodic dimension");
}
// scale factors for CONSTANT and VARIABLE walls
for (int m = 0; m < nwall; m++)
int flag = 0; if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->dimension == 2)
for (int m = 0; m < nwall; m++) error->all(FLERR,
if (wallstyle[m] != EDGE) flag = 1; "Cannot use fix wall/reflect zlo/zhi for a 2d simulation");
if (flag) { // scale factors for CONSTANT and VARIABLE walls
if (scaleflag) {
xscale = domain->lattice->xlattice; int flag = 0;
yscale = domain->lattice->ylattice; for (int m = 0; m < nwall; m++)
zscale = domain->lattice->zlattice; if (wallstyle[m] != EDGE) flag = 1;
}
else xscale = yscale = zscale = 1.0; if (flag) {
if (scaleflag) {
for (int m = 0; m < nwall; m++) { xscale = domain->lattice->xlattice;
if (wallstyle[m] != CONSTANT) continue; yscale = domain->lattice->ylattice;
if (wallwhich[m] < YLO) coord0[m] *= xscale; zscale = domain->lattice->zlattice;
else if (wallwhich[m] < ZLO) coord0[m] *= yscale; }
else coord0[m] *= zscale; else xscale = yscale = zscale = 1.0;
}
} for (int m = 0; m < nwall; m++) {
if (wallstyle[m] != CONSTANT) continue;
// set varflag if any wall positions are variable if (wallwhich[m] < YLO) coord0[m] *= xscale;
else if (wallwhich[m] < ZLO) coord0[m] *= yscale;
varflag = 0; else coord0[m] *= zscale;
for (int m = 0; m < nwall; m++) }
if (wallstyle[m] == VARIABLE) varflag = 1; }
}
} // set varflag if any wall positions are variable
/* ---------------------------------------------------------------------- */ varflag = 0;
for (int m = 0; m < nwall; m++)
FixWallReflect::~FixWallReflect() if (wallstyle[m] == VARIABLE) varflag = 1;
{ }
if (copymode) return;
for (int m = 0; m < nwall; m++)
if (wallstyle[m] == VARIABLE) delete [] varstr[m]; /* ---------------------------------------------------------------------- */
}
FixWallReflect::~FixWallReflect()
/* ---------------------------------------------------------------------- */ {
if (copymode) return;
int FixWallReflect::setmask()
{ for (int m = 0; m < nwall; m++)
int mask = 0; if (wallstyle[m] == VARIABLE) delete [] varstr[m];
mask |= POST_INTEGRATE; }
mask |= POST_INTEGRATE_RESPA;
return mask; /* ---------------------------------------------------------------------- */
}
int FixWallReflect::setmask()
/* ---------------------------------------------------------------------- */ {
int mask = 0;
void FixWallReflect::init() mask |= POST_INTEGRATE;
{ mask |= POST_INTEGRATE_RESPA;
for (int m = 0; m < nwall; m++) { return mask;
if (wallstyle[m] != VARIABLE) continue; }
varindex[m] = input->variable->find(varstr[m]);
if (varindex[m] < 0) /* ---------------------------------------------------------------------- */
error->all(FLERR,"Variable name for fix wall/reflect does not exist");
if (!input->variable->equalstyle(varindex[m])) void FixWallReflect::init()
error->all(FLERR,"Variable for fix wall/reflect is invalid style"); {
} for (int m = 0; m < nwall; m++) {
if (wallstyle[m] != VARIABLE) continue;
int nrigid = 0; varindex[m] = input->variable->find(varstr[m]);
for (int i = 0; i < modify->nfix; i++) if (varindex[m] < 0)
if (modify->fix[i]->rigid_flag) nrigid++; error->all(FLERR,"Variable name for fix wall/reflect does not exist");
if (!input->variable->equalstyle(varindex[m]))
if (nrigid && comm->me == 0) error->all(FLERR,"Variable for fix wall/reflect is invalid style");
error->warning(FLERR,"Should not allow rigid bodies to bounce off " }
"relecting walls");
} int nrigid = 0;
for (int i = 0; i < modify->nfix; i++)
/* ---------------------------------------------------------------------- */ if (modify->fix[i]->rigid_flag) nrigid++;
void FixWallReflect::post_integrate() if (nrigid && comm->me == 0)
{ error->warning(FLERR,"Should not allow rigid bodies to bounce off "
//int m; "relecting walls");
double coord; }
// coord = current position of wall /* ---------------------------------------------------------------------- */
// evaluate variable if necessary, wrap with clear/add
void FixWallReflect::post_integrate()
if (varflag) modify->clearstep_compute(); {
//int m;
for (int m = 0; m < nwall; m++) { double coord;
if (wallstyle[m] == VARIABLE) {
coord = input->variable->compute_equal(varindex[m]); // coord = current position of wall
if (wallwhich[m] < YLO) coord *= xscale; // evaluate variable if necessary, wrap with clear/add
else if (wallwhich[m] < ZLO) coord *= yscale;
else coord *= zscale;
} else coord = coord0[m];
if (varflag) modify->clearstep_compute();
wall_particle(m,wallwhich[m],coord);
for (int m = 0; m < nwall; m++) {
if (varflag) modify->addstep_compute(update->ntimestep + 1); if (wallstyle[m] == VARIABLE) {
} coord = input->variable->compute_equal(varindex[m]);
} if (wallwhich[m] < YLO) coord *= xscale;
else if (wallwhich[m] < ZLO) coord *= yscale;
void FixWallReflect::wall_particle(int m, int which, double coord) else coord *= zscale;
{ } else coord = coord0[m];
int i, dim, side,sign;
double **x = atom->x; wall_particle(m,wallwhich[m],coord);
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal; if (varflag) modify->addstep_compute(update->ntimestep + 1);
}
dim = which / 2; }
side = which % 2;
void FixWallReflect::wall_particle(int m, int which, double coord)
for (i = 0; i < nlocal; i++) { {
if (mask[i] & groupbit) { int i, dim, side,sign;
if (side == 0) { double **x = atom->x;
if (x[i][dim] < coord) { double **v = atom->v;
x[i][dim] = coord + (coord - x[i][dim]); int *mask = atom->mask;
v[i][dim] = -v[i][dim]; int nlocal = atom->nlocal;
}
} else { dim = which / 2;
if (x[i][dim] > coord) { side = which % 2;
x[i][dim] = coord - (x[i][dim] - coord);
v[i][dim] = -v[i][dim]; for (i = 0; i < nlocal; i++)
} if (mask[i] & groupbit) {
}
} if (side == 0) {
} if (x[i][dim] < coord) {
} x[i][dim] = coord + (coord - x[i][dim]);
v[i][dim] = -v[i][dim];
}
} else {
if (x[i][dim] > coord) {
x[i][dim] = coord - (x[i][dim] - coord);
v[i][dim] = -v[i][dim];
}
}
}
}