From 28f5ad69939236d0877267cb436f6d93691d53df Mon Sep 17 00:00:00 2001 From: To Quy-Dong <55282601+toquydong@users.noreply.github.com> Date: Fri, 8 Nov 2019 10:53:07 +0100 Subject: [PATCH] Add files via upload --- src/USER-MISC/fix_wall_reflect_stochastic.cpp | 266 ++++++++++++++++++ src/USER-MISC/fix_wall_reflect_stochastic.h | 85 ++++++ 2 files changed, 351 insertions(+) create mode 100644 src/USER-MISC/fix_wall_reflect_stochastic.cpp create mode 100644 src/USER-MISC/fix_wall_reflect_stochastic.h diff --git a/src/USER-MISC/fix_wall_reflect_stochastic.cpp b/src/USER-MISC/fix_wall_reflect_stochastic.cpp new file mode 100644 index 0000000000..dda2a9296f --- /dev/null +++ b/src/USER-MISC/fix_wall_reflect_stochastic.cpp @@ -0,0 +1,266 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "fix_wall_reflect_stochastic.h" +#include +#include +#include "atom.h" +#include "comm.h" +#include "update.h" +#include "modify.h" +#include "domain.h" +#include "lattice.h" +#include "input.h" +#include "variable.h" +#include "error.h" +#include "force.h" +#include "random_mars.h" + + +using namespace LAMMPS_NS; +using namespace FixConst; + +enum{NONE=0,DIFFUSIVE=1,MAXWELL=2,CERCIGNANILAMPIS=3}; + + +/* ---------------------------------------------------------------------- */ + +FixWallReflectStochastic::FixWallReflectStochastic(LAMMPS *lmp, int narg, char **arg) : + FixWallReflect(lmp, narg, arg), random(NULL) +{ + int arginc,dir; + + if (narg < 8) error->all(FLERR,"Illegal fix wall/reflect/stochastic command"); + + if (domain->triclinic != 0) error->all(FLERR, "This fix wall/reflect/stochastic is not incompatible with triclinic simulation box"); + + dynamic_group_allow = 1; + + // parse args + + nwall = 0; + int scaleflag = 1; + reflectionstyle = NONE; + + + if (strcmp(arg[3],"diffusive") == 0) { + reflectionstyle = DIFFUSIVE; + arginc = 6; + } else if (strcmp(arg[3],"maxwell") == 0) { + reflectionstyle = MAXWELL; + arginc = 7; + } else if (strcmp(arg[3],"cercignanilampis") == 0) { + reflectionstyle = CERCIGNANILAMPIS; + arginc = 9; + } else error->all(FLERR,"Illegal fix wall/reflect/stochastic command"); + + + seedfix = force->inumeric(FLERR,arg[4]); + if (seedfix <=0) error->all(FLERR,"Random seed must be a postive number"); + random = new RanMars(lmp,seedfix + comm->me); + + int iarg = 5; + while (iarg < narg) { + if ((strcmp(arg[iarg],"xlo") == 0) || (strcmp(arg[iarg],"xhi") == 0) || + (strcmp(arg[iarg],"ylo") == 0) || (strcmp(arg[iarg],"yhi") == 0) || + (strcmp(arg[iarg],"zlo") == 0) || (strcmp(arg[iarg],"zhi") == 0)) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix wall/reflect/stochastic command"); + + int newwall; + if (strcmp(arg[iarg],"xlo") == 0) newwall = XLO; + else if (strcmp(arg[iarg],"xhi") == 0) newwall = XHI; + else if (strcmp(arg[iarg],"ylo") == 0) newwall = YLO; + else if (strcmp(arg[iarg],"yhi") == 0) newwall = YHI; + else if (strcmp(arg[iarg],"zlo") == 0) newwall = ZLO; + else if (strcmp(arg[iarg],"zhi") == 0) newwall = ZHI; + + for (int m = 0; (m < nwall) && (m < 6); m++) + if (newwall == wallwhich[m]) + error->all(FLERR,"Wall defined twice in fix wall/reflect/stochastic command"); + + wallwhich[nwall] = newwall; + if (strcmp(arg[iarg+1],"EDGE") == 0) { + wallstyle[nwall] = EDGE; + int dim = wallwhich[nwall] / 2; + int side = wallwhich[nwall] % 2; + if (side == 0) coord0[nwall] = domain->boxlo[dim]; + else coord0[nwall] = domain->boxhi[dim]; + } else { + wallstyle[nwall] = CONSTANT; + coord0[nwall] = force->numeric(FLERR,arg[iarg+1]); + } + + + walltemp[nwall]= force->numeric(FLERR,arg[iarg+2]); + + for (dir = 0; dir < 3; dir++) { + wallvel[nwall][dir]= force->numeric(FLERR,arg[iarg+dir+3]); + int dim = wallwhich[nwall] / 2; + if ((wallvel[nwall][dir] !=0) & (dir == dim)) error->all(FLERR,"The wall velocity must be tangential"); + + if (reflectionstyle == CERCIGNANILAMPIS) { + wallaccom[nwall][dir]= force->numeric(FLERR,arg[iarg+dir+6]); + } else if (reflectionstyle == MAXWELL) { + wallaccom[nwall][dir]= force->numeric(FLERR,arg[iarg+6]);// One accommodation coefficient for all directions + } + } + + nwall++; + iarg += arginc; + + } else if (strcmp(arg[iarg],"units") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal wall/reflect/stochastic command"); + if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; + else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; + else error->all(FLERR,"Illegal fix wall/reflect/stochastic command"); + iarg += 2; + } else error->all(FLERR,"Illegal fix wall/reflect/stochastic command"); + } + + // error check + + if (nwall == 0) error->all(FLERR,"Illegal fix wall 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/stochastic in periodic dimension"); + if ((wallwhich[m] == YLO || wallwhich[m] == YHI) && domain->yperiodic) + error->all(FLERR,"Cannot use fix wall/reflect/stochastic in periodic dimension"); + if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->zperiodic) + error->all(FLERR,"Cannot use fix wall/reflect/stochastic in periodic dimension"); + } + + for (int m = 0; m < nwall; m++) + if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->dimension == 2) + error->all(FLERR, + "Cannot use fix wall/reflect/stochastic zlo/zhi for a 2d simulation"); + + // scale factors for CONSTANT walls, VARIABLE wall is not allowed + + int flag = 0; + for (int m = 0; m < nwall; m++) + if (wallstyle[m] != EDGE) flag = 1; + + if (flag) { + if (scaleflag) { + xscale = domain->lattice->xlattice; + yscale = domain->lattice->ylattice; + zscale = domain->lattice->zlattice; + } + else xscale = yscale = zscale = 1.0; + + for (int m = 0; m < nwall; m++) { + if (wallstyle[m] != CONSTANT) continue; + if (wallwhich[m] < YLO) coord0[m] *= xscale; + else if (wallwhich[m] < ZLO) coord0[m] *= yscale; + else coord0[m] *= zscale; + } + } + + } + + +FixWallReflectStochastic::~FixWallReflectStochastic() +{ + +delete random ; + +} + +/* ---------------------------------------------------------------------- */ + +void FixWallReflectStochastic::wall_particle(int m, int which, double coord) +{ + int i, dir, dim, side, sign; + double vsave,factor,timecol,difftest,theta; + + double *rmass; + double *mass = atom->mass; + + // coord = current position of wall + // evaluate variable if necessary, wrap with clear/add + + double **x = atom->x; + double **v = atom->v; + int *mask = atom->mask; + int *type = atom->type; + int nlocal = atom->nlocal; + + dim = which / 2; + side = which % 2; + + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + + sign = 0; + + if ((side == 0) & (x[i][dim] < coord)){ + sign = 1; + } else if ((side == 1) & (x[i][dim] > coord)){ + sign = -1; + } + + + if (sign != 0) { + if (rmass) theta = force->boltz*walltemp[m]/(rmass[i]*force->mvv2e);// theta = kT/m + else theta = force->boltz*walltemp[m]/(mass[type[i]]*force->mvv2e); + factor = sqrt(theta); + + timecol= (x[i][dim] - coord)/v[i][dim]; // time travelling after collision + + if (reflectionstyle == MAXWELL) {difftest = random->uniform(); }// for Maxwell model only + + for (dir = 0; dir < 3; dir++) { + + x[i][dir] -= v[i][dir]*timecol; // pushing back atoms to the wall position and assign new reflected velocity + + // Diffusive reflection + if (reflectionstyle == DIFFUSIVE){ + + if (dir != dim) { + v[i][dir] = wallvel[m][dir] + random->gaussian(0,factor); + } else { v[i][dir] = sign*random->rayleigh(factor) ; } + + } + + // Maxwell reflection + if (reflectionstyle == MAXWELL){ + + if (difftest < wallaccom[m][dir]) { + if (dir != dim) { + v[i][dir] = wallvel[m][dir] + random->gaussian(0,factor); + } else { v[i][dir] = sign*random->rayleigh(factor) ; } + } else { + if (dir == dim) {v[i][dir] = -v[i][dir]; + } + } + } + + + // Cercignani Lampis reflection + if (reflectionstyle == CERCIGNANILAMPIS){ + if (dir != dim) { + v[i][dir] = wallvel[m][dir] + random->gaussian((1-wallaccom[m][dir])*(v[i][dir]-wallvel[m][dir]),sqrt((2-wallaccom[m][dir])*wallaccom[m][dir]*theta)) ; + } else { v[i][dir] = random->besselexp(theta, wallaccom[m][dir], v[i][dir]) ; } + + } + + // update new position due to the new velocity + if (dir != dim) { + x[i][dir] += v[i][dir]*timecol; } + else {x[i][dir] = coord + v[i][dir]*timecol;} + }// for loop + }// if sign + }// if mask + } +} \ No newline at end of file diff --git a/src/USER-MISC/fix_wall_reflect_stochastic.h b/src/USER-MISC/fix_wall_reflect_stochastic.h new file mode 100644 index 0000000000..c6957bb533 --- /dev/null +++ b/src/USER-MISC/fix_wall_reflect_stochastic.h @@ -0,0 +1,85 @@ +/* -*- c++ -*- ---------------------------------------------------------- + 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 FIX_CLASS + +FixStyle(wall/reflect/stochastic,FixWallReflectStochastic) + +#else + +#ifndef LMP_FIX_WALL_REFLECT_STOCHASTIC_H +#define LMP_FIX_WALL_REFLECT_STOCHASTIC_H + + +#include "random_mars.h" +#include "fix_wall_reflect.h" + +namespace LAMMPS_NS { + +class FixWallReflectStochastic : public FixWallReflect { + public: + FixWallReflectStochastic(class LAMMPS *, int, char **); + void wall_particle(int m,int which, double coord); + virtual ~FixWallReflectStochastic(); + + + + + + protected: + class RanMars *random; + int seedfix; + double walltemp[6],wallvel[6][3],wallaccom[6][3]; + int reflectionstyle; + + +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Wall defined twice in fix wall/stochastic command + +Self-explanatory. + +E: Cannot use fix wall/stochastic in periodic dimension + +Self-explanatory. + +E: Cannot use fix wall/stochastic zlo/zhi for a 2d simulation + +Self-explanatory. + +E: Variable name for fix wall/stochastic does not exist + +Self-explanatory. + +E: Variable for fix wall/stochastic is invalid style + +Only equal-style variables can be used. + +W: Should not allow rigid bodies to bounce off relecting walls + +LAMMPS allows this, but their dynamics are not computed correctly. + +*/