forked from lijiext/lammps
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@10993 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
parent
c3ac2c5be0
commit
83534f0802
|
|
@ -18,8 +18,8 @@ PACKAGE = asphere body class2 colloid dipole fld gpu granular kim \
|
|||
rigid shock srd voronoi xtc
|
||||
|
||||
PACKUSER = user-atc user-awpmd user-cg-cmm user-colvars \
|
||||
user-cuda user-eff user-misc user-omp user-molfile user-phonon \
|
||||
user-reaxc user-sph
|
||||
user-cuda user-eff user-lb user-misc user-omp user-molfile \
|
||||
user-phonon user-reaxc user-sph
|
||||
|
||||
PACKLIB = gpu kim meam poems reax voronoi \
|
||||
user-atc user-awpmd user-colvars user-cuda user-molfile
|
||||
|
|
|
|||
|
|
@ -0,0 +1,54 @@
|
|||
This package contains a LAMMPS implementation of a background
|
||||
Lattice-Boltzmann fluid, which can be used to model MD particles
|
||||
influenced by hydrodynamic forces. Details are described in these two
|
||||
papers:
|
||||
|
||||
Mackay, F. E., Ollila, S.T.T., and Denniston, C., Hydrodynamic Forces
|
||||
Implemented into LAMMPS through a lattice-Boltzmann fluid, Computer
|
||||
Physics Communications 184 (2013) 2021-2031.
|
||||
|
||||
Mackay, F. E., and Denniston, C., Coupling MD particles to a
|
||||
lattice-Boltzmann fluid through the use of conservative forces,
|
||||
J. Comput. Phys. 237 (2013) 289-298.
|
||||
|
||||
See the doc page for the fix lb/fluid command to get started, and see
|
||||
brief descriptions of other fixes below, each of which have their own
|
||||
doc page.
|
||||
|
||||
There are example scripts for using this package in
|
||||
examples/USER/lb.
|
||||
|
||||
IMPORTANT NOTE: This package can only be used if LAMMPS is compiled
|
||||
with MPI (i.e. the serial makefile should not be used to compile the
|
||||
code). Also, several of the test examples provided make use of the
|
||||
rigid fix. Therefore, this should be included in the LAMMPS build by
|
||||
typing "make yes-rigid" prior to the usual compilation (see the
|
||||
"Including/excluding packages" section of the LAMMPS manual).
|
||||
|
||||
The creators of this package are as follows:
|
||||
|
||||
Frances Mackay
|
||||
University of Western Ontario
|
||||
fmackay@uwo.ca
|
||||
|
||||
Dr. Colin Denniston
|
||||
University of Western Ontario
|
||||
cdennist@uwo.ca
|
||||
|
||||
--------------------------------------------------------------------------
|
||||
|
||||
Fixes provided by this package:
|
||||
|
||||
fix_lb_fluid.cpp: fix used to create the lattice-Boltzmann fluid on a
|
||||
grid covering the LAMMPS simulation domain.
|
||||
|
||||
fix_momentum_lb.cpp: fix used to subtract off the total (atom plus fluid)
|
||||
linear momentum from the system.
|
||||
|
||||
fix_pc.cpp: integration algorithm for individual atoms.
|
||||
|
||||
fix_rigid_pc_sphere.cpp: integration algorithm for rigid spherical
|
||||
collections of atoms.
|
||||
|
||||
fix_viscous_lb.cpp: fix to add the fluid force to the atoms when using a
|
||||
built-in LAMMPS integrator.
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -0,0 +1,172 @@
|
|||
/* -----------------------------------------------------------------------
|
||||
LAMMPS 2003 (July 31) - Molecular Dynamics Simulator
|
||||
Sandia National Laboratories, www.cs.sandia.gov/~sjplimp/lammps.html
|
||||
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.
|
||||
|
||||
For more info, see the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------ */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing authors: Frances Mackay, Santtu Ollila, Colin Denniston (UWO)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef FIX_CLASS
|
||||
|
||||
FixStyle(lb/fluid,FixLbFluid)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef FIX_LB_FLUID_H
|
||||
#define FIX_LB_FLUID_H
|
||||
|
||||
#include "stdio.h"
|
||||
#include "fix.h"
|
||||
#include "random_park.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixLbFluid : public Fix {
|
||||
friend class FixLbMomentum;
|
||||
friend class FixLbRigidPCSphere;
|
||||
friend class FixLbPC;
|
||||
friend class FixLbViscous;
|
||||
|
||||
public:
|
||||
FixLbFluid(class LAMMPS *, int, char **);
|
||||
~FixLbFluid();
|
||||
int setmask();
|
||||
void init();
|
||||
void initial_integrate(int);
|
||||
void setup(int);
|
||||
void post_force(int);
|
||||
void end_of_step();
|
||||
|
||||
void grow_arrays(int);
|
||||
void copy_arrays(int, int, int);
|
||||
int pack_exchange(int, double *);
|
||||
int unpack_exchange(int, double *);
|
||||
|
||||
private:
|
||||
#define kappa_lb 0.0
|
||||
#define sqrt2 1.41421356237310
|
||||
|
||||
double viscosity,densityinit_real,a_0_real,T;
|
||||
int setdx,seta0;
|
||||
int numvel;
|
||||
|
||||
double dm_lb,dx_lb,dt_lb; // Lattice units for mass, distance, time.
|
||||
|
||||
int Nbx,Nby,Nbz; // Total # of x,y,z grid points.
|
||||
int subNbx,subNby,subNbz; // # of x,y,z, grid points (including buffer)
|
||||
// on local processor.
|
||||
int me, nprocs; // MPI variables: processor ID, # of processors
|
||||
MPI_Datatype oneslice; // MPI datatypes to pass arrays.
|
||||
MPI_Datatype passxu,passyu,passzu;
|
||||
MPI_Datatype passxf,passyf,passzf;
|
||||
MPI_Datatype passxrho,passyrho,passzrho;
|
||||
MPI_Datatype passxtemp,passytemp,passztemp;
|
||||
|
||||
double kB,densityinit,a_0; // Boltzmann constant, initial density,
|
||||
// and a_0 all in lattice units.
|
||||
double *Gamma;
|
||||
double *NodeArea;
|
||||
int setGamma,setArea;
|
||||
double **hydroF;
|
||||
|
||||
int groupbit_viscouslb, groupbit_pc, groupbit_rigid_pc_sphere;
|
||||
|
||||
double ***density_lb; // fluid density
|
||||
double ****u_lb; // fluid velocity
|
||||
double ****f_lb; // distributions
|
||||
double ****fnew; // used in the calculation of the new
|
||||
// distributions.
|
||||
double ****feq; // equilibrium distributions
|
||||
double ****feqold; // equilibrium distributions from previous
|
||||
// timestep
|
||||
|
||||
double ****feqn; // equilibrium distributions without noise.
|
||||
double ****feqoldn; // equilibrium distributions from previous
|
||||
// timestep without noise.
|
||||
double ****Ff; // Force from the MD particles on the fluid.
|
||||
double ****Fftempx;
|
||||
double ****Fftempy;
|
||||
double ****Fftempz;
|
||||
|
||||
double *Ng_lb; // Lattice Boltzmann variables.
|
||||
double *w_lb;
|
||||
double **mg_lb;
|
||||
int **e;
|
||||
double tau;
|
||||
double expminusdtovertau;
|
||||
double Dcoeff;
|
||||
double K_0;
|
||||
double dtoverdtcollision;
|
||||
|
||||
int step;
|
||||
|
||||
double ****buf; // arrays used to output data.
|
||||
double ****buf2;
|
||||
double ****altogether;
|
||||
double ****altogether2;
|
||||
|
||||
double bodyforcex,bodyforcey,bodyforcez; // Body Forces acting on the fluid (default=0)
|
||||
double vwtp,vwbt; // Velocities of the z walls in the y
|
||||
// direction. (must have fixed boundary
|
||||
// conditions in z)
|
||||
|
||||
int noisestress; // 1 to include noise in the system,
|
||||
// 0 otherwise.
|
||||
double namp,noisefactor;
|
||||
int seed;
|
||||
class RanMars *random;
|
||||
|
||||
int force_diagnostic; // 1 to print out the force action on a group
|
||||
// of particles, 0 otherwise.
|
||||
int igroupforce; // the group for which the force is to be
|
||||
// printed.
|
||||
|
||||
int typeLB;
|
||||
|
||||
int trilinear_stencil; // 1 to use the trilinear stencil, 0 to use the
|
||||
// peskin stencil.
|
||||
|
||||
int readrestart; // 1 to read in data from a restart file.
|
||||
MPI_File pFileRead;
|
||||
|
||||
int printrestart; // 1 to write data to a restart file.
|
||||
MPI_File pFileWrite;
|
||||
|
||||
int printfluid;
|
||||
int fixviscouslb;
|
||||
|
||||
void rescale(void);
|
||||
void (FixLbFluid::*initializeLB)(void);
|
||||
void initializeLB15(void);
|
||||
void initializeLB19(void);
|
||||
void initialize_feq(void);
|
||||
void (FixLbFluid::*equilibriumdist)(int,int,int,int,int,int);
|
||||
void equilibriumdist15(int,int,int,int,int,int);
|
||||
void equilibriumdist19(int,int,int,int,int,int);
|
||||
void parametercalc_part(int,int,int,int,int,int);
|
||||
void parametercalc_full(void);
|
||||
void update_periodic(int,int,int,int,int,int);
|
||||
void (FixLbFluid::*update_full)(void);
|
||||
void update_full15(void);
|
||||
void update_full19(void);
|
||||
void streamout(void);
|
||||
void read_restartfile(void);
|
||||
void write_restartfile(void);
|
||||
void peskin_interpolation(int);
|
||||
void trilinear_interpolation(int);
|
||||
void calc_fluidforce(void);
|
||||
|
||||
};
|
||||
}
|
||||
#endif
|
||||
#endif
|
||||
|
||||
|
|
@ -0,0 +1,296 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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: Frances Mackay, Santtu Ollila, Colin Denniston (UWO)
|
||||
|
||||
Based on fix_momentum,
|
||||
Contributing author: Naveen Michaud-Agrawal (Johns Hopkins U)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "stdlib.h"
|
||||
#include "string.h"
|
||||
#include "fix_lb_momentum.h"
|
||||
#include "atom.h"
|
||||
#include "domain.h"
|
||||
#include "group.h"
|
||||
#include "error.h"
|
||||
#include "fix_lb_fluid.h"
|
||||
#include "modify.h"
|
||||
#include "comm.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixLbMomentum::FixLbMomentum(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg)
|
||||
{
|
||||
if (narg < 4) error->all(FLERR,"Illegal fix lb/momentum command");
|
||||
nevery = atoi(arg[3]);
|
||||
if (nevery <= 0) error->all(FLERR,"Illegal fix lb/momentum command");
|
||||
|
||||
linear = 1;
|
||||
xflag = 1;
|
||||
yflag = 1;
|
||||
zflag = 1;
|
||||
|
||||
int iarg = 4;
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"linear") == 0) {
|
||||
if (iarg+4 > narg) error->all(FLERR,"Illegal fix lb/momentum command");
|
||||
linear = 1;
|
||||
xflag = atoi(arg[iarg+1]);
|
||||
yflag = atoi(arg[iarg+2]);
|
||||
zflag = atoi(arg[iarg+3]);
|
||||
iarg += 4;
|
||||
} else error->all(FLERR,"Illegal fix lb/momentum command");
|
||||
}
|
||||
|
||||
if (linear == 0)
|
||||
error->all(FLERR,"Illegal fix lb/momentum command");
|
||||
|
||||
if (linear)
|
||||
if (xflag < 0 || xflag > 1 || yflag < 0 || yflag > 1 ||
|
||||
zflag < 0 || zflag > 1) error->all(FLERR,"Illegal fix lb/momentum command");
|
||||
|
||||
// cannot have 0 atoms in group
|
||||
|
||||
if (group->count(igroup) == 0.0)
|
||||
error->all(FLERR,"Fix lb/momentum group has no atoms");
|
||||
|
||||
for(int ifix=0; ifix<modify->nfix; ifix++)
|
||||
if(strcmp(modify->fix[ifix]->style,"lb/fluid")==0)
|
||||
fix_lb_fluid = (FixLbFluid *)modify->fix[ifix];
|
||||
|
||||
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixLbMomentum::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= END_OF_STEP;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixLbMomentum::init()
|
||||
{
|
||||
masstotal = group->mass(igroup);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixLbMomentum::end_of_step()
|
||||
{
|
||||
double masslb,masslbloc;
|
||||
double momentumlbloc[3],momentumlb[3];
|
||||
double vcmtotal[3];
|
||||
int numvel = fix_lb_fluid->numvel;
|
||||
double etacov[numvel];
|
||||
double rho;
|
||||
|
||||
if (linear) {
|
||||
double vcm[3];
|
||||
group->vcm(igroup,masstotal,vcm);
|
||||
|
||||
// adjust velocities by vcm to zero linear momentum
|
||||
// only adjust a component if flag is set
|
||||
|
||||
double **v = atom->v;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
int typeLB = fix_lb_fluid->typeLB;
|
||||
int subNbx = fix_lb_fluid->subNbx;
|
||||
int subNby = fix_lb_fluid->subNby;
|
||||
int subNbz = fix_lb_fluid->subNbz;
|
||||
double ***density_lb = fix_lb_fluid->density_lb;
|
||||
double ****f_lb = fix_lb_fluid->f_lb;
|
||||
double ****u_lb = fix_lb_fluid->u_lb;
|
||||
double dm_lb = fix_lb_fluid->dm_lb;
|
||||
double dx_lb = fix_lb_fluid->dx_lb;
|
||||
double dt_lb = fix_lb_fluid->dt_lb;
|
||||
double *Ng_lb = fix_lb_fluid->Ng_lb;
|
||||
double *w_lb = fix_lb_fluid->w_lb;
|
||||
double **mg_lb = fix_lb_fluid->mg_lb;
|
||||
double ucmx,ucmy,ucmz;
|
||||
|
||||
//Calculate the total fluid mass and momentum.
|
||||
masslbloc = 0.0;
|
||||
momentumlbloc[0] = momentumlbloc[1] = momentumlbloc[2] = 0.0;
|
||||
|
||||
for(int i = 1; i<subNbx-1; i++)
|
||||
for(int j = 1; j<subNby-1; j++)
|
||||
for(int k = 1; k<subNbz-1; k++){
|
||||
masslbloc += density_lb[i][j][k];
|
||||
momentumlbloc[0] += density_lb[i][j][k]*u_lb[i][j][k][0];
|
||||
momentumlbloc[1] += density_lb[i][j][k]*u_lb[i][j][k][1];
|
||||
momentumlbloc[2] += density_lb[i][j][k]*u_lb[i][j][k][2];
|
||||
}
|
||||
MPI_Allreduce(&masslbloc,&masslb,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
MPI_Allreduce(&momentumlbloc[0],&momentumlb[0],3,MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
momentumlb[0] *= dm_lb*dx_lb/dt_lb;
|
||||
momentumlb[1] *= dm_lb*dx_lb/dt_lb;
|
||||
momentumlb[2] *= dm_lb*dx_lb/dt_lb;
|
||||
masslb *= dm_lb;
|
||||
|
||||
//Calculate the center of mass velocity of the entire system.
|
||||
vcmtotal[0] = (masstotal*vcm[0] + momentumlb[0])/(masslb + masstotal);
|
||||
vcmtotal[1] = (masstotal*vcm[1] + momentumlb[1])/(masslb + masstotal);
|
||||
vcmtotal[2] = (masstotal*vcm[2] + momentumlb[2])/(masslb + masstotal);
|
||||
|
||||
//Subtract vcm from the particles.
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
if (xflag) v[i][0] -= vcmtotal[0];
|
||||
if (yflag) v[i][1] -= vcmtotal[1];
|
||||
if (zflag) v[i][2] -= vcmtotal[2];
|
||||
}
|
||||
|
||||
vcmtotal[0] *= dt_lb/dx_lb;
|
||||
vcmtotal[1] *= dt_lb/dx_lb;
|
||||
vcmtotal[2] *= dt_lb/dx_lb;
|
||||
|
||||
ucmx = ucmy = ucmz = 0.0;
|
||||
double density_old;
|
||||
double u_old[3];
|
||||
int **e = fix_lb_fluid->e;
|
||||
|
||||
//Subtract vcm from the fluid.
|
||||
for(int i=0; i<subNbx; i++)
|
||||
for(int j=0; j<subNby; j++)
|
||||
for(int k=0; k<subNbz; k++){
|
||||
rho = density_lb[i][j][k];
|
||||
if(xflag) ucmx = vcmtotal[0];
|
||||
if(yflag) ucmy = vcmtotal[1];
|
||||
if(zflag) ucmz = vcmtotal[2];
|
||||
if(numvel==15){
|
||||
etacov[0]=0.0;
|
||||
etacov[1]=rho*ucmx;
|
||||
etacov[2]=rho*ucmy;
|
||||
etacov[3]=rho*ucmz;
|
||||
etacov[4]=rho*(2.*u_lb[i][j][k][0]*ucmx-ucmx*ucmx);
|
||||
etacov[5]=rho*(2.*u_lb[i][j][k][1]*ucmy-ucmy*ucmy);
|
||||
etacov[6]=rho*(2.*u_lb[i][j][k][2]*ucmz-ucmz*ucmz);
|
||||
etacov[7]=rho*(u_lb[i][j][k][0]*ucmy+u_lb[i][j][k][1]*ucmx-ucmx*ucmy);
|
||||
etacov[8]=rho*(u_lb[i][j][k][1]*ucmz+u_lb[i][j][k][2]*ucmy-ucmy*ucmz);
|
||||
etacov[9]=rho*(u_lb[i][j][k][0]*ucmz+u_lb[i][j][k][2]*ucmx-ucmx*ucmz);
|
||||
etacov[10]=0.0;
|
||||
etacov[11]=0.0;
|
||||
etacov[12]=0.0;
|
||||
etacov[13]=rho*(u_lb[i][j][k][0]*u_lb[i][j][k][1]*ucmz+u_lb[i][j][k][0]*ucmy*u_lb[i][j][k][2]-
|
||||
u_lb[i][j][k][0]*ucmy*ucmz+ucmx*u_lb[i][j][k][1]*u_lb[i][j][k][2]-
|
||||
ucmx*u_lb[i][j][k][1]*ucmz-ucmx*ucmy*u_lb[i][j][k][2]+
|
||||
ucmx*ucmy*ucmz);
|
||||
etacov[14]=0.0;
|
||||
}else{
|
||||
etacov[0] = 0.0;
|
||||
etacov[1] = rho*ucmx;
|
||||
etacov[2] = rho*ucmy;
|
||||
etacov[3] = rho*ucmz;
|
||||
etacov[4]=rho*(2.*u_lb[i][j][k][0]*ucmx-ucmx*ucmx);
|
||||
etacov[5]=rho*(2.*u_lb[i][j][k][1]*ucmy-ucmy*ucmy);
|
||||
etacov[6]=rho*(2.*u_lb[i][j][k][2]*ucmz-ucmz*ucmz);
|
||||
etacov[7]=rho*(u_lb[i][j][k][0]*ucmy+u_lb[i][j][k][1]*ucmx-ucmx*ucmy);
|
||||
etacov[8]=rho*(u_lb[i][j][k][0]*ucmz+u_lb[i][j][k][2]*ucmx-ucmx*ucmz);
|
||||
etacov[9]=rho*(u_lb[i][j][k][1]*ucmz+u_lb[i][j][k][2]*ucmy-ucmy*ucmz);
|
||||
etacov[10] = 0.0;
|
||||
etacov[11] = 0.0;
|
||||
etacov[12] = 0.0;
|
||||
etacov[13] = 0.0;
|
||||
etacov[14] = 0.0;
|
||||
etacov[15] = 0.0;
|
||||
etacov[16] = 0.0;
|
||||
etacov[17] = 0.0;
|
||||
etacov[18] = 0.0;
|
||||
}
|
||||
|
||||
for(int l=0; l<numvel; l++)
|
||||
for(int ii=0; ii<numvel; ii++){
|
||||
f_lb[i][j][k][l] -= w_lb[l]*mg_lb[ii][l]*etacov[ii]*Ng_lb[ii];
|
||||
}
|
||||
|
||||
|
||||
if(typeLB == 2){
|
||||
double ****feqold = fix_lb_fluid->feqold;
|
||||
double ****feqoldn = fix_lb_fluid->feqoldn;
|
||||
density_old = 0.0;
|
||||
u_old[0] = u_old[1] = u_old[2] = 0.0;
|
||||
for(int l=0; l<numvel; l++){
|
||||
density_old += feqold[i][j][k][l];
|
||||
u_old[0] += feqold[i][j][k][l]*e[l][0];
|
||||
u_old[1] += feqold[i][j][k][l]*e[l][1];
|
||||
u_old[2] += feqold[i][j][k][l]*e[l][2];
|
||||
}
|
||||
|
||||
if(numvel==15){
|
||||
etacov[0]=0.0;
|
||||
etacov[1]=density_old*ucmx;
|
||||
etacov[2]=density_old*ucmy;
|
||||
etacov[3]=density_old*ucmz;
|
||||
etacov[4]=density_old*(2.*u_old[0]*ucmx-ucmx*ucmx);
|
||||
etacov[5]=density_old*(2.*u_old[1]*ucmy-ucmy*ucmy);
|
||||
etacov[6]=density_old*(2.*u_old[2]*ucmz-ucmz*ucmz);
|
||||
etacov[7]=density_old*(u_old[0]*ucmy+u_old[1]*ucmx-ucmx*ucmy);
|
||||
etacov[8]=density_old*(u_old[1]*ucmz+u_old[2]*ucmy-ucmy*ucmz);
|
||||
etacov[9]=density_old*(u_old[0]*ucmz+u_old[2]*ucmx-ucmx*ucmz);
|
||||
etacov[10]=0.0;
|
||||
etacov[11]=0.0;
|
||||
etacov[12]=0.0;
|
||||
etacov[13]=density_old*(u_old[0]*u_old[1]*ucmz+u_old[0]*ucmy*u_old[2]-
|
||||
u_old[0]*ucmy*ucmz+ucmx*u_old[1]*u_old[2]-
|
||||
ucmx*u_old[1]*ucmz-ucmx*ucmy*u_old[2]+
|
||||
ucmx*ucmy*ucmz);
|
||||
etacov[14]=0.0;
|
||||
}else{
|
||||
etacov[0] = 0.0;
|
||||
etacov[1] = density_old*ucmx;
|
||||
etacov[2] = density_old*ucmy;
|
||||
etacov[3] = density_old*ucmz;
|
||||
etacov[4] = density_old*(2.*u_lb[i][j][k][0]*ucmx-ucmx*ucmx);
|
||||
etacov[5] = density_old*(2.*u_lb[i][j][k][1]*ucmy-ucmy*ucmy);
|
||||
etacov[6] = density_old*(2.*u_lb[i][j][k][2]*ucmz-ucmz*ucmz);
|
||||
etacov[7] = density_old*(u_lb[i][j][k][0]*ucmy+u_lb[i][j][k][1]*ucmx-ucmx*ucmy);
|
||||
etacov[8] = density_old*(u_lb[i][j][k][0]*ucmz+u_lb[i][j][k][2]*ucmx-ucmx*ucmz);
|
||||
etacov[9] = density_old*(u_lb[i][j][k][1]*ucmz+u_lb[i][j][k][2]*ucmy-ucmy*ucmz);
|
||||
etacov[10] = 0.0;
|
||||
etacov[11] = 0.0;
|
||||
etacov[12] = 0.0;
|
||||
etacov[13] = 0.0;
|
||||
etacov[14] = 0.0;
|
||||
etacov[15] = 0.0;
|
||||
etacov[16] = 0.0;
|
||||
etacov[17] = 0.0;
|
||||
etacov[18] = 0.0;
|
||||
}
|
||||
|
||||
for(int l=0; l<numvel; l++)
|
||||
for(int ii=0; ii<numvel; ii++){
|
||||
feqold[i][j][k][l] -= w_lb[l]*mg_lb[ii][l]*etacov[ii]*Ng_lb[ii];
|
||||
feqoldn[i][j][k][l] -= w_lb[l]*mg_lb[ii][l]*etacov[ii]*Ng_lb[ii];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
fix_lb_fluid->parametercalc_full();
|
||||
|
||||
}
|
||||
|
|
@ -0,0 +1,52 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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: Frances Mackay, Santtu Ollila, Colin Denniston (UWO)
|
||||
|
||||
Based on fix_momentum,
|
||||
Contributing author: Naveen Michaud-Agrawal (Johns Hopkins U)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef FIX_CLASS
|
||||
|
||||
FixStyle(lb/momentum,FixLbMomentum)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef FIX_LB_MOMENTUM_H
|
||||
#define FIX_LB_MOMENTUM_H
|
||||
|
||||
#include "fix.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixLbMomentum : public Fix {
|
||||
public:
|
||||
FixLbMomentum(class LAMMPS *, int, char **);
|
||||
int setmask();
|
||||
void init();
|
||||
void end_of_step();
|
||||
|
||||
private:
|
||||
int linear;
|
||||
int xflag,yflag,zflag;
|
||||
double masstotal;
|
||||
|
||||
class FixLbFluid *fix_lb_fluid;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
|
@ -0,0 +1,472 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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: Frances Mackay, Santtu Ollila, Colin Denniston (UWO)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "math.h"
|
||||
#include "stdio.h"
|
||||
#include "string.h"
|
||||
#include "fix_lb_pc.h"
|
||||
#include "atom.h"
|
||||
#include "force.h"
|
||||
#include "update.h"
|
||||
#include "respa.h"
|
||||
#include "error.h"
|
||||
#include "memory.h"
|
||||
#include "comm.h"
|
||||
#include "domain.h"
|
||||
#include "fix_lb_fluid.h"
|
||||
#include "modify.h"
|
||||
#include "mpi.h"
|
||||
#include "group.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixLbPC::FixLbPC(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg)
|
||||
{
|
||||
if (narg < 3)
|
||||
error->all(FLERR,"Illegal fix lb/pc command");
|
||||
|
||||
time_integrate = 1;
|
||||
|
||||
// perform initial allocation of atom-based array
|
||||
// register with Atom class
|
||||
|
||||
force_old = NULL;
|
||||
up = NULL;
|
||||
up_old = NULL;
|
||||
grow_arrays(atom->nmax);
|
||||
atom->add_callback(0);
|
||||
|
||||
Gamma_MD = new double[atom->ntypes+1];
|
||||
|
||||
int groupbit_lb_fluid = 0;
|
||||
for(int ifix=0; ifix<modify->nfix; ifix++)
|
||||
if(strcmp(modify->fix[ifix]->style,"lb/fluid")==0){
|
||||
fix_lb_fluid = (FixLbFluid *)modify->fix[ifix];
|
||||
groupbit_lb_fluid = group->bitmask[modify->fix[ifix]->igroup];
|
||||
}
|
||||
|
||||
if(groupbit_lb_fluid == 0)
|
||||
error->all(FLERR,"the lb/fluid fix must also be used if using the lb/pc fix");
|
||||
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
int print_warning = 0;
|
||||
for(int j=0; j<nlocal; j++){
|
||||
if((mask[j] & groupbit) && !(mask[j] & groupbit_lb_fluid))
|
||||
error->one(FLERR,"can only use the lb/pc fix for an atom if also using the lb/fluid fix for that atom");
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixLbPC::~FixLbPC() {
|
||||
|
||||
atom->delete_callback(id,0);
|
||||
|
||||
memory->destroy(force_old);
|
||||
memory->destroy(up);
|
||||
memory->destroy(up_old);
|
||||
|
||||
delete [] Gamma_MD;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixLbPC::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= INITIAL_INTEGRATE;
|
||||
mask |= FINAL_INTEGRATE;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixLbPC::init()
|
||||
{
|
||||
double *Gamma = fix_lb_fluid->Gamma;
|
||||
double dm_lb = fix_lb_fluid->dm_lb;
|
||||
double dt_lb = fix_lb_fluid->dt_lb;
|
||||
|
||||
MPI_Comm_rank(world,&me);
|
||||
|
||||
dtv = update->dt;
|
||||
dtf = update->dt * force->ftm2v;
|
||||
|
||||
for(int i=0; i<=atom->ntypes; i++)
|
||||
Gamma_MD[i] = Gamma[i]*dm_lb/dt_lb;
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void FixLbPC::initial_integrate(int vflag) {
|
||||
|
||||
double dtfm;
|
||||
|
||||
double **x = atom->x;
|
||||
double dx[3];
|
||||
double **v = atom->v;
|
||||
double **f = atom->f;
|
||||
double *mass = atom->mass;
|
||||
double *rmass = atom->rmass;
|
||||
int *type = atom->type;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
|
||||
|
||||
compute_up();
|
||||
|
||||
for(int i=0; i<nlocal; i++){
|
||||
up_old[i][0] = up[i][0];
|
||||
up_old[i][1] = up[i][1];
|
||||
up_old[i][2] = up[i][2];
|
||||
force_old[i][0] = f[i][0];
|
||||
force_old[i][1] = f[i][1];
|
||||
force_old[i][2] = f[i][2];
|
||||
}
|
||||
|
||||
|
||||
if(rmass){
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
dtfm = dtf/rmass[i];
|
||||
expminusdttimesgamma = exp(-dtv*Gamma_MD[type[i]]/rmass[i]);
|
||||
|
||||
dx[0] = dtv*v[i][0] + 0.5*(f[i][0]*force->ftm2v - Gamma_MD[type[i]]*(v[i][0]-up[i][0]))*dtv*dtv/rmass[i];
|
||||
dx[1] = dtv*v[i][1] + 0.5*(f[i][1]*force->ftm2v - Gamma_MD[type[i]]*(v[i][1]-up[i][1]))*dtv*dtv/rmass[i];
|
||||
dx[2] = dtv*v[i][2] + 0.5*(f[i][2]*force->ftm2v - Gamma_MD[type[i]]*(v[i][2]-up[i][2]))*dtv*dtv/rmass[i];
|
||||
|
||||
x[i][0] += dx[0];
|
||||
x[i][1] += dx[1];
|
||||
x[i][2] += dx[2];
|
||||
|
||||
// Approximation for v
|
||||
if(Gamma_MD[type[i]] == 0.0){
|
||||
v[i][0] += f[i][0]*dtfm;
|
||||
v[i][1] += f[i][1]*dtfm;
|
||||
v[i][2] += f[i][2]*dtfm;
|
||||
}else{
|
||||
v[i][0] = (v[i][0]-up[i][0]-f[i][0]*force->ftm2v/Gamma_MD[type[i]])*expminusdttimesgamma +
|
||||
f[i][0]*force->ftm2v/Gamma_MD[type[i]] + up[i][0];
|
||||
v[i][1] = (v[i][1]-up[i][1]-f[i][1]*force->ftm2v/Gamma_MD[type[i]])*expminusdttimesgamma +
|
||||
f[i][1]*force->ftm2v/Gamma_MD[type[i]] + up[i][1];
|
||||
v[i][2] = (v[i][2]-up[i][2]-f[i][2]*force->ftm2v/Gamma_MD[type[i]])*expminusdttimesgamma +
|
||||
f[i][2]*force->ftm2v/Gamma_MD[type[i]] + up[i][2];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
} else {
|
||||
// this does NOT take varying masses into account
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
dtfm = dtf/mass[type[i]];
|
||||
expminusdttimesgamma = exp(-dtv*Gamma_MD[type[i]]/mass[type[i]]);
|
||||
|
||||
dx[0] = dtv*v[i][0] + 0.5*(f[i][0]*force->ftm2v - Gamma_MD[type[i]]*(v[i][0]-up[i][0]))*dtv*dtv/mass[type[i]];
|
||||
dx[1] = dtv*v[i][1] + 0.5*(f[i][1]*force->ftm2v - Gamma_MD[type[i]]*(v[i][1]-up[i][1]))*dtv*dtv/mass[type[i]];
|
||||
dx[2] = dtv*v[i][2] + 0.5*(f[i][2]*force->ftm2v - Gamma_MD[type[i]]*(v[i][2]-up[i][2]))*dtv*dtv/mass[type[i]];
|
||||
|
||||
x[i][0] += dx[0];
|
||||
x[i][1] += dx[1];
|
||||
x[i][2] += dx[2];
|
||||
|
||||
// Approximation for v
|
||||
if(Gamma_MD[type[i]] == 0.0){
|
||||
v[i][0] += f[i][0]*dtfm;
|
||||
v[i][1] += f[i][1]*dtfm;
|
||||
v[i][2] += f[i][2]*dtfm;
|
||||
}else{
|
||||
v[i][0] = (v[i][0]-up[i][0]-f[i][0]*force->ftm2v/Gamma_MD[type[i]])*expminusdttimesgamma +
|
||||
f[i][0]*force->ftm2v/Gamma_MD[type[i]] + up[i][0];
|
||||
v[i][1] = (v[i][1]-up[i][1]-f[i][1]*force->ftm2v/Gamma_MD[type[i]])*expminusdttimesgamma +
|
||||
f[i][1]*force->ftm2v/Gamma_MD[type[i]] + up[i][1];
|
||||
v[i][2] = (v[i][2]-up[i][2]-f[i][2]*force->ftm2v/Gamma_MD[type[i]])*expminusdttimesgamma +
|
||||
f[i][2]*force->ftm2v/Gamma_MD[type[i]] + up[i][2];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void FixLbPC::final_integrate()
|
||||
{
|
||||
double dtfm;
|
||||
double **v = atom->v;
|
||||
double **f = atom->f;
|
||||
double *mass = atom->mass;
|
||||
double *rmass = atom->rmass;
|
||||
int *type = atom->type;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
|
||||
|
||||
compute_up();
|
||||
|
||||
if(rmass){
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
dtfm = dtf/rmass[i];
|
||||
expminusdttimesgamma = exp(-dtv*Gamma_MD[type[i]]/rmass[i]);
|
||||
DMDcoeff = (dtv - rmass[i]*(1.0-expminusdttimesgamma)/Gamma_MD[type[i]]);
|
||||
|
||||
if(Gamma_MD[type[i]] == 0.0){
|
||||
v[i][0] += 0.5*(f[i][0] - force_old[i][0])*dtfm;
|
||||
v[i][1] += 0.5*(f[i][1] - force_old[i][1])*dtfm;
|
||||
v[i][2] += 0.5*(f[i][2] - force_old[i][2])*dtfm;
|
||||
}else{
|
||||
v[i][0] += DMDcoeff*((f[i][0] - force_old[i][0])*force->ftm2v/Gamma_MD[type[i]] + up[i][0] - up_old[i][0])/dtv;
|
||||
v[i][1] += DMDcoeff*((f[i][1] - force_old[i][1])*force->ftm2v/Gamma_MD[type[i]] + up[i][1] - up_old[i][1])/dtv;
|
||||
v[i][2] += DMDcoeff*((f[i][2] - force_old[i][2])*force->ftm2v/Gamma_MD[type[i]] + up[i][2] - up_old[i][2])/dtv;
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
} else {
|
||||
// this does NOT take varying masses into account
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
dtfm = dtf/mass[type[i]];
|
||||
expminusdttimesgamma = exp(-dtv*Gamma_MD[type[i]]/mass[type[i]]);
|
||||
DMDcoeff = (dtv - mass[type[i]]*(1.0-expminusdttimesgamma)/Gamma_MD[type[i]]);
|
||||
|
||||
if(Gamma_MD[type[i]] == 0.0){
|
||||
v[i][0] += 0.5*(f[i][0] - force_old[i][0])*dtfm;
|
||||
v[i][1] += 0.5*(f[i][1] - force_old[i][1])*dtfm;
|
||||
v[i][2] += 0.5*(f[i][2] - force_old[i][2])*dtfm;
|
||||
}else{
|
||||
v[i][0] += DMDcoeff*((f[i][0] - force_old[i][0])*force->ftm2v/Gamma_MD[type[i]] + up[i][0] - up_old[i][0])/dtv;
|
||||
v[i][1] += DMDcoeff*((f[i][1] - force_old[i][1])*force->ftm2v/Gamma_MD[type[i]] + up[i][1] - up_old[i][1])/dtv;
|
||||
v[i][2] += DMDcoeff*((f[i][2] - force_old[i][2])*force->ftm2v/Gamma_MD[type[i]] + up[i][2] - up_old[i][2])/dtv;
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
allocate atom-based array
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixLbPC::grow_arrays(int nmax)
|
||||
{
|
||||
|
||||
memory->grow(force_old,nmax,3,"FixLbPC:force_old");
|
||||
memory->grow(up_old,nmax,3,"FixLbPC:up_old");
|
||||
memory->grow(up,nmax,3,"FixLbPC:up");
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
copy values within local atom-based array
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixLbPC::copy_arrays(int i, int j, int delflag)
|
||||
{
|
||||
|
||||
force_old[j][0] = force_old[i][0];
|
||||
force_old[j][1] = force_old[i][1];
|
||||
force_old[j][2] = force_old[i][2];
|
||||
up_old[j][0] = up_old[i][0];
|
||||
up_old[j][1] = up_old[i][1];
|
||||
up_old[j][2] = up_old[i][2];
|
||||
up[j][0] = up[i][0];
|
||||
up[j][1] = up[i][1];
|
||||
up[j][2] = up[i][2];
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
pack values in local atom-based array for exchange with another proc
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
int FixLbPC::pack_exchange(int i, double *buf)
|
||||
{
|
||||
|
||||
buf[0] = force_old[i][0];
|
||||
buf[1] = force_old[i][1];
|
||||
buf[2] = force_old[i][2];
|
||||
buf[3] = up_old[i][0];
|
||||
buf[4] = up_old[i][1];
|
||||
buf[5] = up_old[i][2];
|
||||
buf[6] = up[i][0];
|
||||
buf[7] = up[i][1];
|
||||
buf[8] = up[i][2];
|
||||
|
||||
return 9;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
unpack values in local atom-based array from exchange with another proc
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
int FixLbPC::unpack_exchange(int nlocal, double *buf)
|
||||
{
|
||||
|
||||
force_old[nlocal][0] = buf[0];
|
||||
force_old[nlocal][1] = buf[1];
|
||||
force_old[nlocal][2] = buf[2];
|
||||
up_old[nlocal][0] = buf[3];
|
||||
up_old[nlocal][1] = buf[4];
|
||||
up_old[nlocal][2] = buf[5];
|
||||
up[nlocal][0] = buf[6];
|
||||
up[nlocal][1] = buf[7];
|
||||
up[nlocal][2] = buf[8];
|
||||
|
||||
return 9;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void FixLbPC::compute_up(void)
|
||||
{
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
double **x = atom->x;
|
||||
int i,k;
|
||||
int ix,iy,iz;
|
||||
int ixp,iyp,izp;
|
||||
double dx1,dy1,dz1;
|
||||
int isten,ii,jj,kk;
|
||||
double r,rsq,weightx,weighty,weightz;
|
||||
double ****u_lb = fix_lb_fluid->u_lb;
|
||||
int subNbx = fix_lb_fluid->subNbx;
|
||||
int subNby = fix_lb_fluid->subNby;
|
||||
int subNbz = fix_lb_fluid->subNbz;
|
||||
double dx_lb = fix_lb_fluid->dx_lb;
|
||||
double dt_lb = fix_lb_fluid->dt_lb;
|
||||
double FfP[64];
|
||||
int trilinear_stencil = fix_lb_fluid->trilinear_stencil;
|
||||
|
||||
for(i=0; i<nlocal; i++){
|
||||
if(mask[i] & groupbit){
|
||||
|
||||
//Calculate nearest leftmost grid point.
|
||||
//Since array indices from 1 to subNb-2 correspond to the
|
||||
// local subprocessor domain (not indices from 0), use the
|
||||
// ceiling value.
|
||||
ix = (int)ceil((x[i][0]-domain->sublo[0])/dx_lb);
|
||||
iy = (int)ceil((x[i][1]-domain->sublo[1])/dx_lb);
|
||||
iz = (int)ceil((x[i][2]-domain->sublo[2])/dx_lb);
|
||||
|
||||
//Calculate distances to the nearest points.
|
||||
dx1 = x[i][0] - (domain->sublo[0] + (ix-1)*dx_lb);
|
||||
dy1 = x[i][1] - (domain->sublo[1] + (iy-1)*dx_lb);
|
||||
dz1 = x[i][2] - (domain->sublo[2] + (iz-1)*dx_lb);
|
||||
|
||||
// Need to convert these to lattice units:
|
||||
dx1 = dx1/dx_lb;
|
||||
dy1 = dy1/dx_lb;
|
||||
dz1 = dz1/dx_lb;
|
||||
|
||||
up[i][0]=0.0; up[i][1]=0.0; up[i][2]=0.0;
|
||||
if(trilinear_stencil==0){
|
||||
isten=0;
|
||||
for(ii=-1; ii<3; ii++){
|
||||
rsq=(-dx1+ii)*(-dx1+ii);
|
||||
|
||||
if(rsq>=4)
|
||||
weightx=0.0;
|
||||
else{
|
||||
r=sqrt(rsq);
|
||||
if(rsq>1){
|
||||
weightx=(5.0-2.0*r-sqrt(-7.0+12.0*r-4.0*rsq))/8.;
|
||||
} else{
|
||||
weightx=(3.0-2.0*r+sqrt(1.0+4.0*r-4.0*rsq))/8.;
|
||||
}
|
||||
}
|
||||
for(jj=-1; jj<3; jj++){
|
||||
rsq=(-dy1+jj)*(-dy1+jj);
|
||||
if(rsq>=4)
|
||||
weighty=0.0;
|
||||
else{
|
||||
r=sqrt(rsq);
|
||||
if(rsq>1){
|
||||
weighty=(5.0-2.0*r-sqrt(-7.0+12.0*r-4.0*rsq))/8.;
|
||||
} else{
|
||||
weighty=(3.0-2.0*r+sqrt(1.0+4.0*r-4.0*rsq))/8.;
|
||||
}
|
||||
}
|
||||
for(kk=-1; kk<3; kk++){
|
||||
rsq=(-dz1+kk)*(-dz1+kk);
|
||||
if(rsq>=4)
|
||||
weightz=0.0;
|
||||
else{
|
||||
r=sqrt(rsq);
|
||||
if(rsq>1){
|
||||
weightz=(5.0-2.0*r-sqrt(-7.0+12.0*r-4.0*rsq))/8.;
|
||||
} else{
|
||||
weightz=(3.0-2.0*r+sqrt(1.0+4.0*r-4.0*rsq))/8.;
|
||||
}
|
||||
}
|
||||
ixp = ix+ii;
|
||||
iyp = iy+jj;
|
||||
izp = iz+kk;
|
||||
|
||||
|
||||
if(ixp==-1) ixp=subNbx+2;
|
||||
if(iyp==-1) iyp=subNby+2;
|
||||
if(izp==-1) izp=subNbz+2;
|
||||
|
||||
FfP[isten] = weightx*weighty*weightz;
|
||||
// interpolated velocity based on delta function.
|
||||
for(k=0; k<3; k++){
|
||||
up[i][k] += u_lb[ixp][iyp][izp][k]*FfP[isten];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}else{
|
||||
FfP[0] = (1.-dx1)*(1.-dy1)*(1.-dz1);
|
||||
FfP[1] = (1.-dx1)*(1.-dy1)*dz1;
|
||||
FfP[2] = (1.-dx1)*dy1*(1.-dz1);
|
||||
FfP[3] = (1.-dx1)*dy1*dz1;
|
||||
FfP[4] = dx1*(1.-dy1)*(1.-dz1);
|
||||
FfP[5] = dx1*(1.-dy1)*dz1;
|
||||
FfP[6] = dx1*dy1*(1.-dz1);
|
||||
FfP[7] = dx1*dy1*dz1;
|
||||
|
||||
ixp = (ix+1);
|
||||
iyp = (iy+1);
|
||||
izp = (iz+1);
|
||||
|
||||
for (k=0; k<3; k++) { // tri-linearly interpolated velocity at node
|
||||
up[i][k] = u_lb[ix][iy][iz][k]*FfP[0]
|
||||
+ u_lb[ix][iy][izp][k]*FfP[1]
|
||||
+ u_lb[ix][iyp][iz][k]*FfP[2]
|
||||
+ u_lb[ix][iyp][izp][k]*FfP[3]
|
||||
+ u_lb[ixp][iy][iz][k]*FfP[4]
|
||||
+ u_lb[ixp][iy][izp][k]*FfP[5]
|
||||
+ u_lb[ixp][iyp][iz][k]*FfP[6]
|
||||
+ u_lb[ixp][iyp][izp][k]*FfP[7];
|
||||
}
|
||||
}
|
||||
for(k=0; k<3; k++)
|
||||
up[i][k] = up[i][k]*dx_lb/dt_lb;
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,65 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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: Frances Mackay, Santtu Ollila, Colin Denniston (UWO)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef FIX_CLASS
|
||||
|
||||
FixStyle(lb/pc,FixLbPC)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef FIX_LB_PC_H
|
||||
#define FIX_LB_PC_H
|
||||
|
||||
#include "fix.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixLbPC : public Fix {
|
||||
public:
|
||||
FixLbPC(class LAMMPS *, int, char **);
|
||||
~FixLbPC();
|
||||
int setmask();
|
||||
void init();
|
||||
void initial_integrate(int);
|
||||
void final_integrate();
|
||||
|
||||
void grow_arrays(int);
|
||||
void copy_arrays(int, int, int);
|
||||
// void set_arrays(int);
|
||||
int pack_exchange(int, double *);
|
||||
int unpack_exchange(int, double *);
|
||||
|
||||
|
||||
private:
|
||||
double dtv,dtf;
|
||||
int me;
|
||||
double *Gamma_MD;
|
||||
double expminusdttimesgamma;
|
||||
double DMDcoeff;
|
||||
|
||||
double **force_old;
|
||||
double **up;
|
||||
double **up_old;
|
||||
|
||||
void compute_up(void);
|
||||
class FixLbFluid *fix_lb_fluid;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -0,0 +1,116 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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: Frances Mackay, Santtu Ollila, Colin Denniston (UWO)
|
||||
|
||||
Based on fix_rigid.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef FIX_CLASS
|
||||
|
||||
FixStyle(lb/rigid/pc/sphere,FixLbRigidPCSphere)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef FIX_LB_RIGID_PC_SPHERE_H
|
||||
#define FIX_LB_RIGID_PC_SPHERE_H
|
||||
|
||||
#include "fix.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixLbRigidPCSphere : public Fix {
|
||||
public:
|
||||
FixLbRigidPCSphere(class LAMMPS *, int, char **);
|
||||
virtual ~FixLbRigidPCSphere();
|
||||
virtual int setmask();
|
||||
virtual void init();
|
||||
virtual void setup(int);
|
||||
virtual void initial_integrate(int);
|
||||
virtual void final_integrate();
|
||||
virtual double compute_scalar();
|
||||
|
||||
double memory_usage();
|
||||
void grow_arrays(int);
|
||||
void copy_arrays(int, int, int);
|
||||
void set_arrays(int);
|
||||
int pack_exchange(int, double *);
|
||||
int unpack_exchange(int, double *);
|
||||
|
||||
void pre_neighbor();
|
||||
int dof(int);
|
||||
void reset_dt();
|
||||
double compute_array(int, int);
|
||||
|
||||
|
||||
private:
|
||||
double **up;
|
||||
double **up_old;
|
||||
double *Gamma_MD;
|
||||
double expminusdttimesgamma,DMDcoeff;
|
||||
double expminusdttimesgammadiv2;
|
||||
double force_factor,torque_factor;
|
||||
|
||||
double dtv,dtf;
|
||||
|
||||
int nbody; // # of rigid bodies
|
||||
int *nrigid; // # of atoms in each rigid body
|
||||
int *nrigid_shell;
|
||||
double *masstotal; // total mass of each rigid body
|
||||
double *masstotal_shell;
|
||||
double *sphereradius;
|
||||
double **xcm; // coords of center-of-mass of each rigid body
|
||||
double **xcm_old;
|
||||
double **vcm; // velocity of center-of-mass of each
|
||||
double **ucm;
|
||||
double **ucm_old;
|
||||
double **fcm; // force on center-of-mass of each
|
||||
double **fcm_old;
|
||||
double **fcm_fluid;
|
||||
double **omega; // angular momentum of each in space coords
|
||||
double **torque; // torque on each rigid body in space coords
|
||||
double **torque_old;
|
||||
double **torque_fluid;
|
||||
double **torque_fluid_old;
|
||||
double **rotate;
|
||||
int *imagebody; // image flags of xcm of each rigid body
|
||||
double **fflag; // flag for on/off of center-of-mass force
|
||||
double **tflag; // flag for on/off of center-of-mass torque
|
||||
|
||||
int *body; // which body each atom is part of (-1 if none)
|
||||
|
||||
double **sum,**all; // work vectors for each rigid body
|
||||
int **remapflag; // PBC remap flags for each rigid body
|
||||
|
||||
double tfactor; // scale factor on temperature of rigid bodies
|
||||
|
||||
int inner_nodes; // ==1 if certain particle are inside the rigid
|
||||
// body and should not interact with the fluid.
|
||||
// ==0 otherwise.
|
||||
int igroupinner; // specifies the particles which are inside the
|
||||
// spherical rigid body, and do not interact with
|
||||
// the fluid.
|
||||
|
||||
void set_xv();
|
||||
void set_v();
|
||||
|
||||
void compute_up();
|
||||
|
||||
class FixLbFluid *fix_lb_fluid;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
|
@ -0,0 +1,146 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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: Frances Mackay, Santtu Ollila, Colin Denniston (UWO)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "math.h"
|
||||
#include "stdlib.h"
|
||||
#include "string.h"
|
||||
#include "fix_lb_viscous.h"
|
||||
#include "atom.h"
|
||||
#include "update.h"
|
||||
#include "respa.h"
|
||||
#include "error.h"
|
||||
#include "fix_lb_fluid.h"
|
||||
#include "modify.h"
|
||||
#include "group.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixLbViscous::FixLbViscous(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg)
|
||||
{
|
||||
if (narg < 3) error->all(FLERR,"Illegal fix lb/viscous command");
|
||||
|
||||
int groupbit_lb_fluid = 0;
|
||||
|
||||
for(int ifix=0; ifix<modify->nfix; ifix++)
|
||||
if(strcmp(modify->fix[ifix]->style,"lb/fluid")==0){
|
||||
fix_lb_fluid = (FixLbFluid *)modify->fix[ifix];
|
||||
groupbit_lb_fluid = group->bitmask[modify->fix[ifix]->igroup];
|
||||
}
|
||||
|
||||
if(groupbit_lb_fluid == 0)
|
||||
error->all(FLERR,"the lb/fluid fix must also be used if using the lb/viscous fix");
|
||||
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
for(int j=0; j<nlocal; j++){
|
||||
if((mask[j] & groupbit) && !(mask[j] & groupbit_lb_fluid))
|
||||
error->one(FLERR,"to apply a fluid force onto an atom, the lb/fluid fix must be used for that atom.");
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixLbViscous::~FixLbViscous()
|
||||
{
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixLbViscous::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= POST_FORCE;
|
||||
mask |= POST_FORCE_RESPA;
|
||||
mask |= MIN_POST_FORCE;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixLbViscous::init()
|
||||
{
|
||||
|
||||
if (strcmp(update->integrate_style,"respa") == 0)
|
||||
nlevels_respa = ((Respa *) update->integrate)->nlevels;
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixLbViscous::setup(int vflag)
|
||||
{
|
||||
if (strcmp(update->integrate_style,"verlet") == 0)
|
||||
post_force(vflag);
|
||||
else {
|
||||
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
|
||||
post_force_respa(vflag,nlevels_respa-1,0);
|
||||
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixLbViscous::min_setup(int vflag)
|
||||
{
|
||||
post_force(vflag);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixLbViscous::post_force(int vflag)
|
||||
{
|
||||
// apply drag force to atoms in group
|
||||
// direction is opposed to velocity vector
|
||||
// magnitude depends on atom type
|
||||
|
||||
double **f = atom->f;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
double **hydroF = fix_lb_fluid->hydroF;
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
f[i][0] += hydroF[i][0];
|
||||
f[i][1] += hydroF[i][1];
|
||||
f[i][2] += hydroF[i][2];
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixLbViscous::post_force_respa(int vflag, int ilevel, int iloop)
|
||||
{
|
||||
if (ilevel == nlevels_respa-1) post_force(vflag);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixLbViscous::min_post_force(int vflag)
|
||||
{
|
||||
post_force(vflag);
|
||||
}
|
||||
|
||||
|
|
@ -0,0 +1,53 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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: Frances Mackay, Santtu Ollila, Colin Denniston (UWO)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef FIX_CLASS
|
||||
|
||||
FixStyle(lb/viscous,FixLbViscous)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_FIX_LB_VISCOUS_H
|
||||
#define LMP_FIX_LB_VISCOUS_H
|
||||
|
||||
#include "fix.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixLbViscous : public Fix {
|
||||
public:
|
||||
FixLbViscous(class LAMMPS *, int, char **);
|
||||
~FixLbViscous();
|
||||
int setmask();
|
||||
void init();
|
||||
void setup(int);
|
||||
void min_setup(int);
|
||||
void post_force(int);
|
||||
void post_force_respa(int, int, int);
|
||||
void min_post_force(int);
|
||||
|
||||
protected:
|
||||
int nlevels_respa;
|
||||
|
||||
private:
|
||||
class FixLbFluid *fix_lb_fluid;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
|
@ -686,7 +686,7 @@ void Special::angle_trim()
|
|||
int nbuf = 0;
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
if (num_angle && atom->nangles) nbuf += 2*num_angle[i];
|
||||
if (num_dihedral && atom->ndihedrals) nbuf + 2*2*num_dihedral[i];
|
||||
if (num_dihedral && atom->ndihedrals) nbuf += 2*2*num_dihedral[i];
|
||||
}
|
||||
int *buf;
|
||||
memory->create(buf,nbuf,"special:buf");
|
||||
|
|
|
|||
|
|
@ -245,6 +245,48 @@ void Update::set_units(const char *style)
|
|||
dt = 0.001;
|
||||
neighbor->skin = 2.0;
|
||||
|
||||
} else if (strcmp(style,"micro") == 0) {
|
||||
force->boltz = 1.3806504e-8;
|
||||
force->hplanck = 6.62606896e-13;
|
||||
force->mvv2e = 1.0;
|
||||
force->ftm2v = 1.0;
|
||||
force->mv2d = 1.0;
|
||||
force->nktv2p = 1.0;
|
||||
force->qqr2e = 8.9876e30;
|
||||
force->qe2f = 1.0;
|
||||
force->vxmu2f = 1.0;
|
||||
force->xxt2kmu = 1.0;
|
||||
force->e_mass = 0.0; // not yet set
|
||||
force->hhmrr2e = 0.0;
|
||||
force->mvh2r = 0.0;
|
||||
force->angstrom = 1.0e-4;
|
||||
force->femtosecond = 1.0e-9;
|
||||
force->qelectron = 1.6021765e-19;
|
||||
|
||||
dt = 2.0;
|
||||
neighbor->skin = 0.1;
|
||||
|
||||
} else if (strcmp(style,"nano") == 0) {
|
||||
force->boltz = 0.013806503;
|
||||
force->hplanck = 6.62606896e-4;
|
||||
force->mvv2e = 1.0;
|
||||
force->ftm2v = 1.0;
|
||||
force->mv2d = 1.0;
|
||||
force->nktv2p = 1.0;
|
||||
force->qqr2e = 8.9876e39;
|
||||
force->qe2f = 1.0;
|
||||
force->vxmu2f = 1.0;
|
||||
force->xxt2kmu = 1.0;
|
||||
force->e_mass = 0.0; // not yet set
|
||||
force->hhmrr2e = 0.0;
|
||||
force->mvh2r = 0.0;
|
||||
force->angstrom = 1.0e-1;
|
||||
force->femtosecond = 1.0e-6;
|
||||
force->qelectron = 1.6021765e-19;
|
||||
|
||||
dt = 0.00045;
|
||||
neighbor->skin = 0.1;
|
||||
|
||||
} else error->all(FLERR,"Illegal units command");
|
||||
|
||||
delete [] unit_style;
|
||||
|
|
|
|||
Loading…
Reference in New Issue