diff --git a/src/Makefile b/src/Makefile index abb32bc41e..0263fbae21 100755 --- a/src/Makefile +++ b/src/Makefile @@ -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 diff --git a/src/USER-LB/README b/src/USER-LB/README new file mode 100755 index 0000000000..290c2f923d --- /dev/null +++ b/src/USER-LB/README @@ -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. diff --git a/src/USER-LB/fix_lb_fluid.cpp b/src/USER-LB/fix_lb_fluid.cpp new file mode 100755 index 0000000000..8b052878de --- /dev/null +++ b/src/USER-LB/fix_lb_fluid.cpp @@ -0,0 +1,3358 @@ +/* ---------------------------------------------------------------------- + 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 "fix_lb_fluid.h" +#include "math.h" +#include "mpi.h" +#include "stdlib.h" +#include "stdio.h" +#include "string.h" +#include "comm.h" +#include "memory.h" +#include "error.h" +#include "domain.h" +#include "atom.h" +#include +#include +#include "group.h" +#include "random_mars.h" +#include "update.h" +#include "force.h" +#include "modify.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +FixLbFluid::FixLbFluid(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + //===================================================================================================== + // Sample inputfile call: + // fix # group lb/fluid nevery typeLB viscosity densityinit_real + // + // where: nevery: call this fix every nevery timesteps. + // (keep this set to 1 for now). + // typeLB: there are two different integrators + // in the code labelled "1" and "2". + // viscosity: the viscosity of the fluid. + // densityinit_real: the density of the fluid. + // + // optional arguments: + // "setArea" type node_area: set the surface area per node associated with a + // given atom type. By default the surface area + // is set at 1.0*dx_lb^2. + // "setGamma" gamma: specify a user-defined value for the force + // coupling constant, instead of using the default + // value. + // "scaleGamma" type scale_factor: scale the user provided force coupling constant + // by the factor, scale_factor, for the given atom + // type. + // "dx" dx_lb: the lattice-Boltzmann grid spacing. + // "dm" dm_lb: the lattice-Boltzmann mass unit. + // "a0" a_0_real: the square of the sound speed in the fluid. + // "noise" Temperature seed: include noise in the system. + // Temperature is the temperature for the fluid. + // seed is the seed for the random number generator. + // "calcforce" N group: print the force acting on a given group every + // N timesteps. + // "trilinear": use the trilinear interpolation stencil. + // "read_restart" restart_file: restart a fluid run from restart_file. + // "write_restart" N: write a fluid restart file every N timesteps. + // "zwall_velocity" velocity_bottom velocity_top: assign velocities to the z-walls + // in the system. + // "bodyforce" bodyforcex bodyforcey bodyforcez: add a constant body force to the + // fluid. + // "printfluid" N: print the fluid density and velocity at each + // grid point every N timesteps. + // "D3Q19": use the 19 velocity D3Q19 model. By default, + // the 15 velocity D3Q15 model is used. + //===================================================================================================== + + if(narg <7) error->all(FLERR,"Illegal fix lb/fluid command"); + + MPI_Comm_rank(world,&me); + MPI_Comm_size(world,&nprocs); + + nevery = atoi(arg[3]); + typeLB = atoi(arg[4]); + viscosity = atof(arg[5]); + densityinit_real = atof(arg[6]); + + // Default values for optional arguments: + force_diagnostic=0; + noisestress = 0; + trilinear_stencil = 0; + readrestart = 0; + printrestart = 0; + bodyforcex = bodyforcey = bodyforcez = 0.0; + vwtp = vwbt = 0.0; + printfluid = 0; + T = 300.0; + dm_lb = 1.0; + fixviscouslb = 0; + setdx = 1; + seta0 = 1; + setGamma = 0; + setArea = 0; + numvel = 15; + + Gamma = NULL; + NodeArea = NULL; + + int iarg = 7; + while (iarg < narg){ + if(strcmp(arg[iarg],"setArea")==0){ + if(setGamma == 1) + error->all(FLERR,"Illegal fix lb/fluid command: cannot use a combination of default and user-specified gamma values"); + setArea = 1; + int itype = atoi(arg[iarg+1]); + double areafactor = atof(arg[iarg+2]); + if(itype <= 0 || itype > atom->ntypes || areafactor < 0.0) + error->all(FLERR,"Illegal fix lb/fluid command: setArea"); + if(NodeArea == NULL){ + NodeArea = new double[atom->ntypes+1]; + for(int i=0; i<=atom->ntypes; i++) NodeArea[i] = -1.0; + } + NodeArea[itype] = areafactor; + iarg += 3; + } + else if(strcmp(arg[iarg],"setGamma")==0){ + if(setArea == 1) + error->all(FLERR,"Illegal fix lb/fluid command: cannot use a combination of default and user-specified gamma values"); + setGamma = 1; + double Gammaone; + Gammaone = atof(arg[iarg+1]); + if(Gamma == NULL) + Gamma = new double[atom->ntypes+1]; + for(int i=0; i<=atom->ntypes; i++) Gamma[i] = Gammaone; + iarg += 2; + } + else if(strcmp(arg[iarg],"scaleGamma")==0){ + if(setGamma == 0) + error->all(FLERR,"Illegal fix lb/fluid command: must set a value for Gamma before scaling it"); + int itype = atoi(arg[iarg+1]); + double scalefactor = atof(arg[iarg+2]); + if(itype <= 0 || itype > atom->ntypes || scalefactor < 0.0) + error->all(FLERR,"Illegal fix lb/fluid command: scaleGamma"); + Gamma[itype] *= scalefactor; + iarg += 3; + } + else if(strcmp(arg[iarg],"dx")==0){ + dx_lb = atof(arg[iarg+1]); + iarg += 2; + setdx = 0; + } + else if(strcmp(arg[iarg],"dm")==0){ + dm_lb = atof(arg[iarg+1]); + iarg += 2; + } + else if(strcmp(arg[iarg],"a0")==0){ + a_0_real = atof(arg[iarg+1]); + iarg += 2; + seta0 = 0; + } + else if(strcmp(arg[iarg],"noise")== 0){ + noisestress = 1; + T = atof(arg[iarg+1]); + seed = atoi(arg[iarg+2]); + iarg += 3; + } + else if(strcmp(arg[iarg],"calcforce")==0){ + force_diagnostic = atoi(arg[iarg+1]); + igroupforce=group->find(arg[iarg+2]); + iarg += 3; + } + else if(strcmp(arg[iarg],"trilinear")==0){ + trilinear_stencil = 1; + iarg += 1; + } + else if(strcmp(arg[iarg],"read_restart")==0){ + readrestart = 1; + int nlength = strlen(arg[iarg+1]) + 16; + char *filename = new char[nlength]; + strcpy(filename,arg[iarg+1]); + MPI_File_open(world,filename,MPI_MODE_RDONLY,MPI_INFO_NULL,&pFileRead); + delete [] filename; + iarg += 2; + } + else if(strcmp(arg[iarg],"write_restart")==0){ + printrestart = atoi(arg[iarg+1]); + iarg += 2; + } + else if(strcmp(arg[iarg],"zwall_velocity")==0){ + if(domain->periodicity[2]!=0) error->all(FLERR,"fix lb/fluid error: setting \ +a z wall velocity without implementing fixed BCs in z"); + vwbt = atof(arg[iarg+1]); + vwtp = atof(arg[iarg+2]); + iarg += 3; + } + else if(strcmp(arg[iarg],"bodyforce")==0){ + bodyforcex = atof(arg[iarg+1]); + bodyforcey = atof(arg[iarg+2]); + bodyforcez = atof(arg[iarg+3]); + iarg += 4; + } + else if(strcmp(arg[iarg],"printfluid")==0){ + printfluid = atoi(arg[iarg+1]); + iarg += 2; + } + else if(strcmp(arg[iarg],"D3Q19")==0){ + numvel = 19; + iarg += 1; + } + else error->all(FLERR,"Illegal fix lb/fluid command"); + } + + //-------------------------------------------------------------------------- + //Choose between D3Q15 and D3Q19 functions: + //-------------------------------------------------------------------------- + if(numvel == 15){ + initializeLB = &FixLbFluid::initializeLB15; + equilibriumdist = &FixLbFluid::equilibriumdist15; + update_full = &FixLbFluid::update_full15; + }else{ + initializeLB = &FixLbFluid::initializeLB19; + equilibriumdist = &FixLbFluid::equilibriumdist19; + update_full = &FixLbFluid::update_full19; + } + + //-------------------------------------------------------------------------- + // perform initial allocation of atom-based array register + // with Atom class + //-------------------------------------------------------------------------- + hydroF = NULL; + grow_arrays(atom->nmax); + atom->add_callback(0); + + for(int i=0; inmax; i++) + for(int j=0; j<3; j++) + hydroF[i][j] = 0.0; + + Ng_lb = NULL; + w_lb = NULL; + mg_lb = NULL; + e = NULL; + feq = NULL; + feqold = NULL; + feqn = NULL; + feqoldn = NULL; + f_lb = NULL; + fnew = NULL; + density_lb = NULL; + u_lb = NULL; + altogether = NULL; + buf = NULL; + Ff = NULL; + Fftempx = NULL; + Fftempy = NULL; + Fftempz = NULL; + + //-------------------------------------------------------------------------- + // Set the lattice Boltzmann dt. + //-------------------------------------------------------------------------- + dt_lb=nevery*(update->dt); + + //-------------------------------------------------------------------------- + // Set the lattice Boltzmann dx if it wasn't specified in the + // input. + //-------------------------------------------------------------------------- + if(setdx == 1){ + double dx_lb1 = sqrt(3.0*viscosity*dt_lb/densityinit_real); + double mindomain = std::min(std::min(domain->xprd/comm->procgrid[0],domain->yprd/comm->procgrid[1]),domain->zprd/comm->procgrid[2]); + dx_lb = mindomain/floor(mindomain/dx_lb1); + + if(comm->me==0){ + char str[128]; + sprintf(str,"Setting the lattice-Boltzmann dx to %10.6f",dx_lb); + error->message(FLERR,str); + } + } + //-------------------------------------------------------------------------- + // If the area per node has not been set by the user, set to the + // default value of dx_lb*dx_lb. + //-------------------------------------------------------------------------- + if(setGamma == 0){ + if(setArea == 0){ + if(comm->me==0){ + error->message(FLERR,"Assuming an area per node of dx*dx for all of the MD particles. This should only be used if these all correspond to point particles; otherwise, change using the setArea keyword"); + } + NodeArea = new double[atom->ntypes+1]; + for(int i=0; i<=atom->ntypes; i++) NodeArea[i] = -1.0; + } + for(int i=0; i<=atom->ntypes; i++) + if(NodeArea[i] < 0.0) NodeArea[i] = dx_lb*dx_lb; + } + //-------------------------------------------------------------------------- + // Set a0 if it wasn't specified in the input + //-------------------------------------------------------------------------- + if(seta0 == 1) + a_0_real = 0.33333333*dx_lb*dx_lb/dt_lb/dt_lb; + + //-------------------------------------------------------------------------- + // Check to make sure that the total number of grid points in each direction + // divides evenly among the processors in that direction. + // Shrink-wrapped boundary conditions (which are not permitted by this fix) + // might cause a problem, so check for this. A full check of the boundary + // conditions is performed in the init routine, rather than here, as it is + // possible to change the BCs between runs. + //-------------------------------------------------------------------------- + double aa; + double eps=1.0e-8; + aa = (domain->xprd/comm->procgrid[0])/dx_lb; + if(fabs(aa - floor(aa+0.5)) > eps){ + if(domain->boundary[0][0] != 0){ + error->all(FLERR,"the x-direction must be periodic"); + } + char errormessage[200]; + sprintf(errormessage,"With dx= %f, and the simulation domain divided by %i processors in the x direction, the simulation domain in the x direction must be a multiple of %f",dx_lb,comm->procgrid[0],comm->procgrid[0]*dx_lb); + error->all(FLERR,errormessage); + } + aa = (domain->yprd/comm->procgrid[1])/dx_lb; + if(fabs(aa - floor(aa+0.5)) > eps){ + if(domain->boundary[1][0] != 0){ + error->all(FLERR,"the y-direction must be periodic"); + } + char errormessage[200]; + sprintf(errormessage,"With dx= %f, and the simulation domain divided by %i processors in the y direction, the simulation domain in the y direction must be a multiple of %f",dx_lb,comm->procgrid[1],comm->procgrid[1]*dx_lb); + error->all(FLERR,errormessage); + } + aa = (domain->zprd/comm->procgrid[2])/dx_lb; + if(fabs(aa - floor(aa+0.5)) > eps){ + if(domain->boundary[2][0] == 2 || domain->boundary[2][0] == 3){ + error->all(FLERR,"the z-direction can not have shrink-wrap boundary conditions"); + } + char errormessage[200]; + sprintf(errormessage,"With dx= %f, and the simulation domain divided by %i processors in the z direction, the simulation domain in the z direction must be a multiple of %f",dx_lb,comm->procgrid[2],comm->procgrid[2]*dx_lb); + error->all(FLERR,errormessage); + } + + //-------------------------------------------------------------------------- + // Set the total number of grid points in each direction. + //-------------------------------------------------------------------------- + Nbx = (int)(domain->xprd/dx_lb + 0.5); + Nby = (int)(domain->yprd/dx_lb + 0.5); + Nbz = (int)(domain->zprd/dx_lb + 0.5); + + //-------------------------------------------------------------------------- + // Set the number of grid points in each dimension for the local subgrids. + //-------------------------------------------------------------------------- + subNbx= Nbx/comm->procgrid[0] + 2; + subNby= Nby/comm->procgrid[1] + 2; + subNbz= Nbz/comm->procgrid[2] + 2; + + //-------------------------------------------------------------------------- + // In order to calculate the fluid forces correctly, need to have atleast + // 5 grid points in each direction per processor. + //-------------------------------------------------------------------------- + if(subNbx<7 || subNby < 7 || subNbz<7) + error->all(FLERR,"Need at least 5 grid points in each direction per processor"); + + // If there are walls in the z-direction add an extra grid point. + if(domain->periodicity[2]==0){ + Nbz += 1; + if(comm->myloc[2]==comm->procgrid[2]-1) + subNbz += 1; + } + + if(comm->me==0){ + char str[128]; + if(setdx == 1){ + sprintf(str,"Using a lattice-Boltzmann grid of %i by %i by %i total grid points. To change, use the dx keyword",Nbx,Nby,Nbz); + }else{ + sprintf(str,"Using a lattice-Boltzmann grid of %i by %i by %i total grid points.",Nbx,Nby,Nbz); + } + error->message(FLERR,str); + } + + //-------------------------------------------------------------------------- + // Store the largest value of subNbz, which is needed for allocating the + // buf array (since a processor with comm->myloc[2] == comm->procgrid[2]-1 + // may have an additional subNbz point as compared with the rest). + //-------------------------------------------------------------------------- + int subNbzmax; + MPI_Allreduce(&subNbz,&subNbzmax,1,MPI_INT,MPI_MAX,world); + + //-------------------------------------------------------------------------- + // Create the MPI datatypes used to pass portions of arrays: + // datatypes to pass the f and feq arrays. + //-------------------------------------------------------------------------- + MPI_Aint sizeofdouble; + MPI_Type_extent(MPI_DOUBLE,&sizeofdouble); + + MPI_Type_vector(subNbz-2,numvel,numvel,MPI_DOUBLE,&oneslice); + MPI_Type_commit(&oneslice); + MPI_Type_hvector(subNby-2,1,numvel*subNbz*sizeofdouble,oneslice,&passxf); + MPI_Type_commit(&passxf); + + MPI_Type_hvector(subNbx,1,numvel*subNbz*subNby*sizeofdouble,oneslice,&passyf); + MPI_Type_commit(&passyf); + + MPI_Type_free(&oneslice); + MPI_Type_vector(subNby,numvel,numvel*subNbz,MPI_DOUBLE,&oneslice); + MPI_Type_commit(&oneslice); + MPI_Type_hvector(subNbx,1,numvel*subNbz*subNby*sizeofdouble,oneslice,&passzf); + MPI_Type_commit(&passzf); + + // datatypes to pass the u array, and the Ff array. + MPI_Type_free(&oneslice); + MPI_Type_vector(subNbz+3,3,3,MPI_DOUBLE,&oneslice); + MPI_Type_commit(&oneslice); + MPI_Type_hvector(subNby+3,1,3*(subNbz+3)*sizeofdouble,oneslice,&passxu); + MPI_Type_commit(&passxu); + + MPI_Type_hvector(subNbx+3,1,3*(subNbz+3)*(subNby+3)*sizeofdouble,oneslice,&passyu); + MPI_Type_commit(&passyu); + + MPI_Type_free(&oneslice); + MPI_Type_vector(subNby+3,3,3*(subNbz+3),MPI_DOUBLE,&oneslice); + MPI_Type_commit(&oneslice); + MPI_Type_hvector(subNbx+3,1,3*(subNbz+3)*(subNby+3)*sizeofdouble,oneslice,&passzu); + MPI_Type_commit(&passzu); + + // datatypes to pass the density array. + MPI_Type_free(&oneslice); + MPI_Type_vector(subNbz+3,1,1,MPI_DOUBLE,&oneslice); + MPI_Type_commit(&oneslice); + MPI_Type_hvector(subNby+3,1,1*(subNbz+3)*sizeofdouble,oneslice,&passxrho); + MPI_Type_commit(&passxrho); + + MPI_Type_hvector(subNbx+3,1,1*(subNbz+3)*(subNby+3)*sizeofdouble,oneslice,&passyrho); + MPI_Type_commit(&passyrho); + + MPI_Type_free(&oneslice); + MPI_Type_vector(subNby+3,1,1*(subNbz+3),MPI_DOUBLE,&oneslice); + MPI_Type_commit(&oneslice); + MPI_Type_hvector(subNbx+3,1,1*(subNbz+3)*(subNby+3)*sizeofdouble,oneslice,&passzrho); + MPI_Type_commit(&passzrho); + + // datatypes to receive a portion of the Ff array. + MPI_Type_free(&oneslice); + MPI_Type_vector(subNbz+3,3,3,MPI_DOUBLE,&oneslice); + MPI_Type_commit(&oneslice); + MPI_Type_hvector(subNby+3,1,3*(subNbz+3)*sizeofdouble,oneslice,&passxtemp); + MPI_Type_commit(&passxtemp); + + MPI_Type_hvector(subNbx+3,1,3*(subNbz+3)*5*sizeofdouble,oneslice,&passytemp); + MPI_Type_commit(&passytemp); + + MPI_Type_free(&oneslice); + MPI_Type_vector(subNby+3,3,3*5,MPI_DOUBLE,&oneslice); + MPI_Type_commit(&oneslice); + MPI_Type_hvector(subNbx+3,1,3*5*(subNby+3)*sizeofdouble,oneslice,&passztemp); + MPI_Type_commit(&passztemp); + + MPI_Type_free(&oneslice); + + //-------------------------------------------------------------------------- + // Allocate the necessary arrays. + //-------------------------------------------------------------------------- + memory->create(Ng_lb,numvel,"FixLbFluid:Ng_lb"); + memory->create(w_lb,numvel,"FixLbFluid:w_lb"); + memory->create(mg_lb,numvel,numvel,"FixLbFluid:mg_lb"); + memory->create(e,numvel,3,"FixLbFluid:e"); + memory->create(feq,subNbx,subNby,subNbz,numvel,"FixLbFluid:feq"); + if(typeLB == 2){ + memory->create(feqold,subNbx,subNby,subNbz,numvel,"FixLbFluid:feqold"); + memory->create(feqn,subNbx,subNby,subNbz,numvel,"FixLbFluid:feqn"); + memory->create(feqoldn,subNbx,subNby,subNbz,numvel,"FixLbFluid:feqoldn"); + } + memory->create(f_lb,subNbx,subNby,subNbz,numvel,"FixLbFluid:f_lb"); + memory->create(fnew,subNbx,subNby,subNbz,numvel,"FixLbFluid:fnew"); + memory->create(density_lb,subNbx+3,subNby+3,subNbz+3,"FixLbFluid:density_lb"); + memory->create(u_lb,subNbx+3,subNby+3,subNbz+3,3,"FixLbFluid:u_lb"); + if(printfluid > 0){ + memory->create(buf,subNbx,subNby,subNbzmax,4,"FixLbFluid:buf"); + if(me==0) + memory->create(altogether,Nbx,Nby,Nbz,4,"FixLbFluid:altogether"); + } + memory->create(Ff,subNbx+3,subNby+3,subNbz+3,3,"FixLbFluid:Ff"); + memory->create(Fftempx,5,subNby+3,subNbz+3,3,"FixLbFluid:Fftempx"); + memory->create(Fftempy,subNbx+3,5,subNbz+3,3,"FixLbFluid:Fftempy"); + memory->create(Fftempz,subNbx+3,subNby+3,5,3,"FixLbFluid:Fftempz"); + + if(noisestress==1){ + random = new RanMars(lmp,seed + comm->me); + } + + //-------------------------------------------------------------------------- + // Rescale the variables to Lattice Boltzmann dimensionless units. + //-------------------------------------------------------------------------- + rescale(); + + //-------------------------------------------------------------------------- + // Initialize the arrays. + //-------------------------------------------------------------------------- + (*this.*initializeLB)(); + initialize_feq(); + +} + +FixLbFluid::~FixLbFluid() +{ + + atom->delete_callback(id,0); + memory->destroy(hydroF); + + memory->destroy(Ng_lb); + memory->destroy(w_lb); + memory->destroy(mg_lb); + memory->destroy(e); + memory->destroy(feq); + if(typeLB == 2){ + memory->destroy(feqold); + memory->destroy(feqn); + memory->destroy(feqoldn); + } + memory->destroy(f_lb); + memory->destroy(fnew); + memory->destroy(density_lb); + memory->destroy(u_lb); + if(printfluid>0){ + if(me==0) + memory->destroy(altogether); + memory->destroy(buf); + } + memory->destroy(Ff); + memory->destroy(Fftempx); + memory->destroy(Fftempy); + memory->destroy(Fftempz); + + if(noisestress==1){ + delete random; + } + + if(setGamma == 1){ + delete [] Gamma; + }else{ + delete [] NodeArea; + } +} + +int FixLbFluid::setmask() +{ + int mask =0; + mask |= INITIAL_INTEGRATE; + mask |= POST_FORCE; + mask |= END_OF_STEP; + return mask; +} + +void FixLbFluid::init(void) +{ + + int i,j; + + //-------------------------------------------------------------------------- + // Check to see if the MD timestep has changed between runs. + //-------------------------------------------------------------------------- + double dt_lb_now; + dt_lb_now=nevery*(update->dt); + + if(fabs(dt_lb_now - dt_lb) > 1.0e-12){ + error->warning(FLERR,"Timestep has changed between runs with the same lb/fluid. Unphysical results may occur"); + } + + //-------------------------------------------------------------------------- + // Make sure the size of the simulation domain has not changed + // between runs. + //-------------------------------------------------------------------------- + int Nbx_now,Nby_now,Nbz_now; + Nbx_now = (int)(domain->xprd/dx_lb + 0.5); + Nby_now = (int)(domain->yprd/dx_lb + 0.5); + Nbz_now = (int)(domain->zprd/dx_lb + 0.5); + // If there are walls in the z-direction add an extra grid point. + if(domain->periodicity[2]==0){ + Nbz_now += 1; + } + + if(Nbx_now != Nbx || Nby_now != Nby || Nbz_now != Nbz){ + error->all(FLERR,"the simulation domain can not change shape between runs with the same lb/fluid"); + } + + //-------------------------------------------------------------------------- + // Check to make sure that the chosen LAMMPS boundary types are compatible + // with this fix. + // shrink-wrap is not compatible in any dimension. + // fixed only works in the z-direction. + //-------------------------------------------------------------------------- + if(domain->boundary[0][0] != 0){ + error->all(FLERR,"the x-direction must be periodic"); + } + if(domain->boundary[1][0] != 0){ + error->all(FLERR,"the y-direction must be periodic"); + } + if(domain->boundary[2][0] == 2 || domain->boundary[2][0] == 3){ + error->all(FLERR,"the z-direction can not have shrink-wrap boundary conditions"); + } + + //-------------------------------------------------------------------------- + // Check if the lb/viscous fix is also called: + //-------------------------------------------------------------------------- + groupbit_viscouslb = groupbit_pc = groupbit_rigid_pc_sphere = 0; + for (i = 0; i < modify->nfix; i++){ + if (strcmp(modify->fix[i]->style,"lb/viscous") == 0){ + fixviscouslb = 1; + groupbit_viscouslb = group->bitmask[modify->fix[i]->igroup]; + } + if(strcmp(modify->fix[i]->style,"lb/pc")==0){ + groupbit_pc = group->bitmask[modify->fix[i]->igroup]; + } + if(strcmp(modify->fix[i]->style,"lb/rigid/pc/sphere")==0){ + groupbit_rigid_pc_sphere = group->bitmask[modify->fix[i]->igroup]; + } + } + + // Warn if the fluid force is not applied to any of the particles. + if(!(groupbit_viscouslb || groupbit_pc || groupbit_rigid_pc_sphere) && comm->me==0){ + error->message(FLERR,"Not adding the fluid force to any of the MD particles. To add this force use one of the lb/viscous, lb/pc, or lb/rigid/pc/sphere fixes"); + } + + // If fix lb/viscous is called for a particular atom, make sure + // lb/pc or lb/rigid/pc/sphere are not: + if(fixviscouslb == 1){ + int *mask = atom->mask; + int nlocal = atom->nlocal; + for(j=0; jone(FLERR,"should not use the lb/viscous command when integrating with the lb/pc fix"); + if((mask[j] & groupbit) && (mask[j] & groupbit_viscouslb) && (mask[j] & groupbit_rigid_pc_sphere)) + error->one(FLERR,"should not use the lb/viscous command when integrating with the lb/rigid/pc/sphere fix"); + } + } + +} + +void FixLbFluid::setup(int vflag) +{ + //-------------------------------------------------------------------------- + // Need to calculate the force on the fluid for a restart run. + //-------------------------------------------------------------------------- + if(step > 0) + calc_fluidforce(); +} + +void FixLbFluid::initial_integrate(int vflag) +{ + //-------------------------------------------------------------------------- + // Print a header labelling any output printed to the screen. + //-------------------------------------------------------------------------- + static int printheader = 1; + + if(printheader == 1){ + if(force_diagnostic > 0 && me == 0){ + printf("-------------------------------------------------------------------------------\n"); + printf(" F_x F_y F_z T_x T_y T_z\n"); + printf("-------------------------------------------------------------------------------\n"); + } + + if(printfluid > 0 && me == 0){ + printf("---------------------------------------------------------------------\n"); + printf(" density u_x u_y u_z \n"); + printf("---------------------------------------------------------------------\n"); + } + printheader = 0; + } + + //-------------------------------------------------------------------------- + // Determine the equilibrium distribution on the local subgrid. + //-------------------------------------------------------------------------- + (*this.*equilibriumdist)(1,subNbx-1,1,subNby-1,1,subNbz-1); + + //-------------------------------------------------------------------------- + // Using the equilibrium distribution, calculate the new + // distribution function. + //-------------------------------------------------------------------------- + (*this.*update_full)(); + + std::swap(f_lb,fnew); + + //-------------------------------------------------------------------------- + // Calculate moments of the distribution function. + //-------------------------------------------------------------------------- + parametercalc_full(); + + //-------------------------------------------------------------------------- + // Store the equilibrium distribution function, it is needed in + // the next time step by the update routine. + //-------------------------------------------------------------------------- + if(typeLB == 2){ + std::swap(feqold,feq); + std::swap(feqoldn,feqn); + } + + //-------------------------------------------------------------------------- + // Perform diagnostics, and print output for the graphics program + //-------------------------------------------------------------------------- + if(printfluid > 0 && update->ntimestep > 0 && (update->ntimestep % printfluid == 0)) + streamout(); + +} +void FixLbFluid::post_force(int vflag) +{ + if(fixviscouslb==1) + calc_fluidforce(); +} + +void FixLbFluid::end_of_step() +{ + if(fixviscouslb==0) + calc_fluidforce(); + + if(printrestart>0){ + if((update->ntimestep)%printrestart == 0){ + write_restartfile(); + } + } + +} + +//========================================================================== +// allocate atom-based array +//========================================================================== +void FixLbFluid::grow_arrays(int nmax) +{ + memory->grow(hydroF,nmax,3,"FixLbFluid:hydroF"); +} + +//========================================================================== +// copy values within local atom-based array +//========================================================================== +void FixLbFluid::copy_arrays(int i, int j, int delflag) +{ + hydroF[j][0] = hydroF[i][0]; + hydroF[j][1] = hydroF[i][1]; + hydroF[j][2] = hydroF[i][2]; +} + +//========================================================================== +// pack values in local atom-based array for exchange with another proc +//========================================================================== +int FixLbFluid::pack_exchange(int i, double *buf) +{ + buf[0] = hydroF[i][0]; + buf[1] = hydroF[i][1]; + buf[2] = hydroF[i][2]; + + return 3; +} + +//========================================================================== +// unpack values in local atom-based array from exchange with another proc +//========================================================================== +int FixLbFluid::unpack_exchange(int nlocal, double *buf) +{ + hydroF[nlocal][0] = buf[0]; + hydroF[nlocal][1] = buf[1]; + hydroF[nlocal][2] = buf[2]; + + return 3; +} + +//========================================================================== +// calculate the force from the local atoms acting on the fluid. +//========================================================================== +void FixLbFluid::calc_fluidforce(void) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + double **x = atom->x; + int i,j,k,m; + MPI_Request requests[20]; + MPI_Status statuses[20]; + double forceloc[3],force[3]; + double torqueloc[3],torque[3]; + int numrequests; + + //-------------------------------------------------------------------------- + // Zero out arrays + //-------------------------------------------------------------------------- + std::fill(&Ff[0][0][0][0],&Ff[0][0][0][0] + (subNbx+3)*(subNby+3)*(subNbz+3)*3,0.0); + std::fill(&Fftempx[0][0][0][0],&Fftempx[0][0][0][0] + 5*(subNby+3)*(subNbz+3)*3,0.0); + std::fill(&Fftempy[0][0][0][0],&Fftempy[0][0][0][0] + (subNbx+3)*5*(subNbz+3)*3,0.0); + std::fill(&Fftempz[0][0][0][0],&Fftempz[0][0][0][0] + (subNbx+3)*(subNby+3)*5*3,0.0); + + forceloc[0] = forceloc[1] = forceloc[2] = 0.0; + torqueloc[0] = torqueloc[1] = torqueloc[2] = 0.0; + + for(i=0; inmax; i++) + for(j=0; j<3; j++) + hydroF[i][j] = 0.0; + + + double unwrap[3]; + double dx,dy,dz; + double massone; + tagint *image = atom->image; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; + double sum[4],xcm[4]; + + if(force_diagnostic > 0 && update->ntimestep > 0 && (update->ntimestep % force_diagnostic == 0)){ + //Calculate the center of mass of the particle group + //(needed to calculate the torque). + sum[0] = sum[1] = sum[2] = sum[3] = 0.0; + for(i=0; ibitmask[igroupforce]){ + + domain->unmap(x[i],image[i],unwrap); + + if(rmass) massone = rmass[i]; + else massone = mass[type[i]]; + + sum[0] += unwrap[0]*massone; + sum[1] += unwrap[1]*massone; + sum[2] += unwrap[2]*massone; + sum[3] += massone; + } + } + MPI_Allreduce(&sum[0],&xcm[0],4,MPI_DOUBLE,MPI_SUM,world); + xcm[0] = xcm[0]/xcm[3]; + xcm[1] = xcm[1]/xcm[3]; + xcm[2] = xcm[2]/xcm[3]; + } + + //-------------------------------------------------------------------------- + //Calculate the contribution to the force on the fluid. + //-------------------------------------------------------------------------- + for(i=0; i 0 && update->ntimestep > 0 && (update->ntimestep % force_diagnostic == 0)){ + if(mask[i] & group->bitmask[igroupforce]){ + + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - xcm[0]; + dy = unwrap[1] - xcm[1]; + dz = unwrap[2] - xcm[2]; + + forceloc[0] += hydroF[i][0]; + forceloc[1] += hydroF[i][1]; + forceloc[2] += hydroF[i][2]; + torqueloc[0] += dy*hydroF[i][2] - dz*hydroF[i][1]; + torqueloc[1] += dz*hydroF[i][0] - dx*hydroF[i][2]; + torqueloc[2] += dx*hydroF[i][1] - dy*hydroF[i][0]; + } + } + } + } + + //-------------------------------------------------------------------------- + //Communicate the force contributions which lie outside the local processor + //sub domain. + //-------------------------------------------------------------------------- + for(i=0; i<10; i++) + requests[i]=MPI_REQUEST_NULL; + MPI_Isend(&Ff[0][0][0][0],1,passxu,comm->procneigh[0][0],10,world,&requests[0]); + MPI_Isend(&Ff[subNbx+2][0][0][0],1,passxu,comm->procneigh[0][0],20,world,&requests[1]); + MPI_Isend(&Ff[subNbx-1][0][0][0],1,passxu,comm->procneigh[0][1],30,world,&requests[2]); + MPI_Isend(&Ff[subNbx][0][0][0],1,passxu,comm->procneigh[0][1],40,world,&requests[3]); + MPI_Isend(&Ff[subNbx+1][0][0][0],1,passxu,comm->procneigh[0][1],50,world,&requests[4]); + MPI_Irecv(&Fftempx[0][0][0][0],1,passxtemp,comm->procneigh[0][1],10,world,&requests[5]); + MPI_Irecv(&Fftempx[1][0][0][0],1,passxtemp,comm->procneigh[0][1],20,world,&requests[6]); + MPI_Irecv(&Fftempx[2][0][0][0],1,passxtemp,comm->procneigh[0][0],30,world,&requests[7]); + MPI_Irecv(&Fftempx[3][0][0][0],1,passxtemp,comm->procneigh[0][0],40,world,&requests[8]); + MPI_Irecv(&Fftempx[4][0][0][0],1,passxtemp,comm->procneigh[0][0],50,world,&requests[9]); + MPI_Waitall(10,requests,statuses); + + for(j=0; jprocneigh[1][0],10,world,&requests[0]); + MPI_Isend(&Ff[0][subNby+2][0][0],1,passyu,comm->procneigh[1][0],20,world,&requests[1]); + MPI_Isend(&Ff[0][subNby-1][0][0],1,passyu,comm->procneigh[1][1],30,world,&requests[2]); + MPI_Isend(&Ff[0][subNby][0][0],1,passyu,comm->procneigh[1][1],40,world,&requests[3]); + MPI_Isend(&Ff[0][subNby+1][0][0],1,passyu,comm->procneigh[1][1],50,world,&requests[4]); + MPI_Irecv(&Fftempy[0][0][0][0],1,passytemp,comm->procneigh[1][1],10,world,&requests[5]); + MPI_Irecv(&Fftempy[0][1][0][0],1,passytemp,comm->procneigh[1][1],20,world,&requests[6]); + MPI_Irecv(&Fftempy[0][2][0][0],1,passytemp,comm->procneigh[1][0],30,world,&requests[7]); + MPI_Irecv(&Fftempy[0][3][0][0],1,passytemp,comm->procneigh[1][0],40,world,&requests[8]); + MPI_Irecv(&Fftempy[0][4][0][0],1,passytemp,comm->procneigh[1][0],50,world,&requests[9]); + MPI_Waitall(10,requests,statuses); + + for(i=0; iprocneigh[2][0],10,world,&requests[0]); + MPI_Isend(&Ff[0][0][subNbz+2][0],1,passzu,comm->procneigh[2][0],20,world,&requests[1]); + MPI_Isend(&Ff[0][0][subNbz-1][0],1,passzu,comm->procneigh[2][1],30,world,&requests[2]); + MPI_Isend(&Ff[0][0][subNbz][0],1,passzu,comm->procneigh[2][1],40,world,&requests[3]); + MPI_Isend(&Ff[0][0][subNbz+1][0],1,passzu,comm->procneigh[2][1],50,world,&requests[4]); + MPI_Irecv(&Fftempz[0][0][0][0],1,passztemp,comm->procneigh[2][1],10,world,&requests[5]); + MPI_Irecv(&Fftempz[0][0][1][0],1,passztemp,comm->procneigh[2][1],20,world,&requests[6]); + MPI_Irecv(&Fftempz[0][0][2][0],1,passztemp,comm->procneigh[2][0],30,world,&requests[7]); + MPI_Irecv(&Fftempz[0][0][3][0],1,passztemp,comm->procneigh[2][0],40,world,&requests[8]); + MPI_Irecv(&Fftempz[0][0][4][0],1,passztemp,comm->procneigh[2][0],50,world,&requests[9]); + MPI_Waitall(10,requests,statuses); + + for(i=0; i 0 && update->ntimestep > 0 && (update->ntimestep % force_diagnostic == 0)){ + force[0] = force[1] = force[2] = 0.0; + torque[0] = torque[1] = torque[2] =0.0; + + MPI_Allreduce(&forceloc[0],&force[0],3,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&torqueloc[0],&torque[0],3,MPI_DOUBLE,MPI_SUM,world); + + if(me==0){ + printf("%E %E %E %E %E %E\n",force[0],force[1],force[2], + torque[0],torque[1],torque[2]); + + } + } + +} +//========================================================================== +// uses the Peskin stencil to perform the velocity, density and +// force interpolations. +//========================================================================== +void FixLbFluid::peskin_interpolation(int i) +{ + double **x = atom->x; + double **v = atom->v; + int *type = atom->type; + double *rmass = atom->rmass; + double *mass = atom->mass; + double massone; + int ix,iy,iz; + int ixp,iyp,izp; + double dx1,dy1,dz1; + int isten,ii,jj,kk; + double r,rsq,weightx,weighty,weightz; + double FfP[64]; + int k; + double unode[3]; + double mnode; + double gammavalue; + + //-------------------------------------------------------------------------- + //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; + + unode[0]=0.0; unode[1]=0.0; unode[2]=0.0; + mnode = 0.0; + isten=0; + + //-------------------------------------------------------------------------- + // Calculate the interpolation weights, and interpolated values of + // the fluid velocity, and density. + //-------------------------------------------------------------------------- + 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; + + //The atom is allowed to be within one lattice grid point outside the + //local processor sub-domain. + if(ixp < -1 || ixp > (subNbx+1) || iyp < -1 || iyp > (subNby+1) || izp < -1 || izp > (subNbz+1)) + error->one(FLERR,"Atom outside local processor simulation domain. Either unstable fluid pararmeters, or \ +require more frequent neighborlist rebuilds"); + + if(domain->periodicity[2] == 0 && comm->myloc[2] == 0 && izp < 1) + error->warning(FLERR,"Atom too close to lower z wall. Unphysical results may occur"); + if(domain->periodicity[2] == 0 && comm->myloc[2] == (comm->procgrid[2]-1) && (izp > (subNbz-2) )) + error->warning(FLERR,"Atom too close to upper z wall. Unphysical results may occur"); + + 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++){ + unode[k] += u_lb[ixp][iyp][izp][k]*FfP[isten]; + } + if(setGamma==0) + mnode += density_lb[ixp][iyp][izp]*FfP[isten]; + + isten++; + } + } + } + if(setGamma==0){ + mnode *= NodeArea[type[i]]; + + if(rmass) massone = rmass[i]; + else massone = mass[type[i]]; + massone = massone/dm_lb; + + gammavalue = 2.0*(mnode*massone)*dtoverdtcollision/(mnode+massone); + } + else{ + gammavalue = Gamma[type[i]]; + } + + isten=0; + for(ii=-1; ii<3; ii++) + for(jj=-1; jj<3; jj++) + for(kk=-1; kk<3; kk++){ + 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; + // Compute the force on the fluid. Need to convert the velocity from + // LAMMPS units to LB units. + for(k=0; k<3; k++){ + Ff[ixp][iyp][izp][k] += gammavalue*((v[i][k]*dt_lb/dx_lb)-unode[k])*FfP[isten]; + } + + isten++; + } + for(k=0; k<3; k++) + hydroF[i][k] = -1.0*gammavalue*((v[i][k]*dt_lb/dx_lb)-unode[k])*dm_lb*dx_lb/dt_lb/dt_lb; +} + +//========================================================================== +// uses the trilinear stencil to perform the velocity, density and +// force interpolations. +//========================================================================== +void FixLbFluid::trilinear_interpolation(int i) +{ + double **x = atom->x; + double **v = atom->v; + int *type = atom->type; + double *rmass = atom->rmass; + double *mass = atom->mass; + double massone; + int ix,iy,iz; + int ixp,iyp,izp; + double dx1,dy1,dz1; + double FfP[8]; + int k; + double unode[3]; + double mnode; + double gammavalue; + + //-------------------------------------------------------------------------- + // 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; + + //-------------------------------------------------------------------------- + // Calculate the interpolation weights + //-------------------------------------------------------------------------- + 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); + + //The atom is allowed to be within one lattice grid point outside the + //local processor sub-domain. + if(ix < 0 || ixp > (subNbx+1) || iy < 0 || iyp > (subNby+1) || iz < 0 || izp > (subNbz+1)) + error->one(FLERR,"Atom outside local processor simulation domain. Either unstable fluid pararmeters, or \ +require more frequent neighborlist rebuilds"); + + if(domain->periodicity[2] == 0 && comm->myloc[2] == 0 && (iz < 1 || izp < 1)) + error->warning(FLERR,"Atom too close to lower z wall. Unphysical results may occur"); + if(domain->periodicity[2] == 0 && comm->myloc[2] == (comm->procgrid[2]-1) && (izp > (subNbz-2) || iz > (subNbz-2))) + error->warning(FLERR,"Atom too close to upper z wall. Unphysical results may occur"); + + + for (k=0; k<3; k++) { // tri-linearly interpolated velocity at node + unode[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]; + } + + if(setGamma==0){ + mnode = density_lb[ix][iy][iz]*FfP[0] + + density_lb[ix][iy][izp]*FfP[1] + + density_lb[ix][iyp][iz]*FfP[2] + + density_lb[ix][iyp][izp]*FfP[3] + + density_lb[ixp][iy][iz]*FfP[4] + + density_lb[ixp][iy][izp]*FfP[5] + + density_lb[ixp][iyp][iz]*FfP[6] + + density_lb[ixp][iyp][izp]*FfP[7]; + + mnode *= NodeArea[type[i]]; + + if(rmass) massone = rmass[i]; + else massone = mass[type[i]]; + massone = massone/dm_lb; + + gammavalue = 2.0*(mnode*massone)*dtoverdtcollision/(mnode+massone); + }else{ + gammavalue = Gamma[type[i]]; + } + + + for(k=0; k<3; k++){ + Ff[ix][iy][iz][k] += gammavalue*((v[i][k]*dt_lb/dx_lb)-unode[k])*FfP[0]; + Ff[ix][iy][izp][k] += gammavalue*((v[i][k]*dt_lb/dx_lb)-unode[k])*FfP[1]; + Ff[ix][iyp][iz][k] += gammavalue*((v[i][k]*dt_lb/dx_lb)-unode[k])*FfP[2]; + Ff[ix][iyp][izp][k] += gammavalue*((v[i][k]*dt_lb/dx_lb)-unode[k])*FfP[3]; + Ff[ixp][iy][iz][k] += gammavalue*((v[i][k]*dt_lb/dx_lb)-unode[k])*FfP[4]; + Ff[ixp][iy][izp][k] += gammavalue*((v[i][k]*dt_lb/dx_lb)-unode[k])*FfP[5]; + Ff[ixp][iyp][iz][k] += gammavalue*((v[i][k]*dt_lb/dx_lb)-unode[k])*FfP[6]; + Ff[ixp][iyp][izp][k] += gammavalue*((v[i][k]*dt_lb/dx_lb)-unode[k])*FfP[7]; + } + + for(k=0; k<3; k++) + hydroF[i][k] = -1.0*gammavalue*((v[i][k]*dt_lb/dx_lb)-unode[k])*dm_lb*dx_lb/dt_lb/dt_lb; + +} + +//========================================================================== +// read in a fluid restart file. This is only used to restart the +// fluid portion of a LAMMPS simulation. +//========================================================================== +void FixLbFluid::read_restartfile(void) +{ + MPI_Status status; + MPI_Datatype realtype; + MPI_Datatype filetype; + + + int realsizes[4] = {subNbx,subNby,subNbz,numvel}; + int realstarts[4] = {1,1,1,0}; + int gsizes[4] = {Nbx,Nby,Nbz,numvel}; + int lsizes[4] = {subNbx-2,subNby-2,subNbz-2,numvel}; + int starts[4] = {comm->myloc[0]*(subNbx-2),comm->myloc[1]*(subNby-2),comm->myloc[2]*(subNbz-2),0}; + if(domain->periodicity[2]==0 && comm->myloc[2]==comm->procgrid[2]-1){ + starts[2] = comm->myloc[2]*(subNbz-3); + } + + MPI_Type_create_subarray(4,realsizes,lsizes,realstarts,MPI_ORDER_C,MPI_DOUBLE,&realtype); + MPI_Type_commit(&realtype); + + MPI_Type_create_subarray(4,gsizes,lsizes,starts,MPI_ORDER_C,MPI_DOUBLE,&filetype); + MPI_Type_commit(&filetype); + + MPI_File_set_view(pFileRead,0,MPI_DOUBLE,filetype,"native",MPI_INFO_NULL); + MPI_File_seek(pFileRead,0,MPI_SEEK_SET); + MPI_File_read_all(pFileRead,&f_lb[0][0][0][0],1,realtype,&status); + if(typeLB == 2){ + MPI_File_read_all(pFileRead,&feqold[0][0][0][0],1,realtype,&status); + MPI_File_read_all(pFileRead,&feqoldn[0][0][0][0],1,realtype,&status); + } + + MPI_Type_free(&realtype); + MPI_Type_free(&filetype); + MPI_File_close(&pFileRead); + +} + +//========================================================================== +// write a fluid restart file. +//========================================================================== +void FixLbFluid::write_restartfile(void) +{ + + MPI_File fh; + MPI_Status status; + MPI_Datatype realtype; + MPI_Datatype filetype; + + char *hfile; + hfile = new char[32]; + sprintf(hfile,"FluidRestart_%d.dat",update->ntimestep); + + MPI_File_open(world,hfile,MPI_MODE_WRONLY | MPI_MODE_CREATE, MPI_INFO_NULL,&fh); + + int realsizes[4] = {subNbx,subNby,subNbz,numvel}; + int realstarts[4] = {1,1,1,0}; + int gsizes[4] = {Nbx,Nby,Nbz,numvel}; + int lsizes[4] = {subNbx-2,subNby-2,subNbz-2,numvel}; + int starts[4] = {comm->myloc[0]*(subNbx-2),comm->myloc[1]*(subNby-2),comm->myloc[2]*(subNbz-2),0}; + if(domain->periodicity[2]==0 && comm->myloc[2]==comm->procgrid[2]-1){ + starts[2] = comm->myloc[2]*(subNbz-3); + } + + MPI_Type_create_subarray(4,realsizes,lsizes,realstarts,MPI_ORDER_C,MPI_DOUBLE,&realtype); + MPI_Type_commit(&realtype); + + MPI_Type_create_subarray(4,gsizes,lsizes,starts,MPI_ORDER_C,MPI_DOUBLE,&filetype); + MPI_Type_commit(&filetype); + + MPI_File_set_view(fh,0,MPI_DOUBLE,filetype,"native",MPI_INFO_NULL); + MPI_File_write_all(fh,&f_lb[0][0][0][0],1,realtype,&status); + if(typeLB == 2){ + MPI_File_write_all(fh,&feqold[0][0][0][0],1,realtype,&status); + MPI_File_write_all(fh,&feqoldn[0][0][0][0],1,realtype,&status); + } + + MPI_Type_free(&realtype); + MPI_Type_free(&filetype); + MPI_File_close(&fh); + delete [] hfile; + +} + +//========================================================================== +// rescale the simulation parameters so that dx_lb=dt_lb=dm_lb=1. +// This assumes that all the simulation parameters have been given in +// terms of distance, time and mass units. +//========================================================================== +void FixLbFluid::rescale(void) +{ + vwtp = vwtp*dt_lb/dx_lb; + vwbt = vwbt*dt_lb/dx_lb; + + bodyforcex = bodyforcex*dt_lb*dt_lb/dx_lb; + bodyforcey = bodyforcey*dt_lb*dt_lb/dx_lb; + bodyforcez = bodyforcez*dt_lb*dt_lb/dx_lb; + + tau=(3.0*viscosity/densityinit_real)*dt_lb*dt_lb/dx_lb/dx_lb; + tau /= dt_lb; + if(typeLB==1) + tau = tau + 0.5; + + if(setGamma == 0){ + for(int i=0; i<= atom->ntypes; i++){ + NodeArea[i] = NodeArea[i]/dx_lb/dx_lb; + } + }else{ + for(int i=0; i<= atom->ntypes; i++){ + Gamma[i] = Gamma[i]*dt_lb/dm_lb; + } + } + + densityinit = densityinit_real*dx_lb*dx_lb*dx_lb/dm_lb; + + a_0 = a_0_real*dt_lb*dt_lb/(dx_lb*dx_lb); + + // Warn if using the D3Q19 model with noise, and a0 is too small. + if(numvel==19 && noisestress==1 && a_0 < 0.2){ + error->warning(FLERR,"Fix lb/fluid WARNING: Chosen value for a0 may be too small. \ +Check temperature reproduction.\n"); + } + + if(noisestress==1){ + if(a_0>0.5555555){ + error->all(FLERR,"Fix lb/fluid ERROR: the Lattice Boltzmann dx and dt need \ +to be chosen such that the scaled a_0 < 5/9\n"); + } + } + + // Courant Condition: + if(a_0 >= 1.0){ + error->all(FLERR,"Fix lb/fluid ERROR: the lattice Boltzmann dx and dt do not \ +satisfy the Courant condition.\n"); + } + + kB = (force->boltz/force->mvv2e)*dt_lb*dt_lb/dx_lb/dx_lb/dm_lb; + + if(typeLB==1){ + expminusdtovertau = 0.0; + Dcoeff = 0.0; + namp = 2.0*kB*T*(tau-0.5)/3.0; + noisefactor = 1.0; + if(a_0 <= 0.333333333333333){ + K_0 = 5.17*(0.333333333333333 - a_0); + }else{ + K_0 = 2.57*(a_0 - 0.333333333333333); + } + dtoverdtcollision = dt_lb*6.0*viscosity/densityinit_real/dx_lb/dx_lb; + }else if(typeLB==2){ + expminusdtovertau=exp(-1.0/tau); + Dcoeff=(1.0-(1.0-expminusdtovertau)*tau); + namp = 2.0*kB*T/3.; + noisefactor=sqrt((1.0-expminusdtovertau*expminusdtovertau)/ + (2.0))/(1.0-expminusdtovertau); + K_0 = 4.5*(1.0/3.0-a_0); + dtoverdtcollision = dt_lb*3.0*viscosity/densityinit_real/dx_lb/dx_lb; + } + +} + +//========================================================================== +// Set the lattice-Boltzmann velocity vectors and weights for the D3Q15 +// model. Initialize the fluid velocity and density. +//========================================================================== +void FixLbFluid::initializeLB15(void) +{ + int i,j,k,m; + + //velocity vectors. + e[0][0]= 0; + e[0][1]= 0; + e[0][2]= 0; + + e[1][0]= 1; + e[1][1]= 0; + e[1][2]= 0; + + e[2][0]= 0; + e[2][1]= 1; + e[2][2]= 0; + + e[3][0]= -1; + e[3][1]= 0; + e[3][2]= 0; + + e[4][0]= 0; + e[4][1]= -1; + e[4][2]= 0; + + e[5][0]= 0; + e[5][1]= 0; + e[5][2]= 1; + + e[6][0]= 0; + e[6][1]= 0; + e[6][2]= -1; + + e[7][0]= 1; + e[7][1]= 1; + e[7][2]= 1; + + e[8][0]= -1; + e[8][1]= 1; + e[8][2]= 1; + + e[9][0]= -1; + e[9][1]= -1; + e[9][2]= 1; + + e[10][0]= 1; + e[10][1]= -1; + e[10][2]= 1; + + e[11][0]= 1; + e[11][1]= 1; + e[11][2]= -1; + + e[12][0]= -1; + e[12][1]= 1; + e[12][2]= -1; + + e[13][0]= -1; + e[13][1]= -1; + e[13][2]= -1; + + e[14][0]= 1; + e[14][1]= -1; + e[14][2]= -1; + + //weights. + w_lb[0]=2./9.; + w_lb[1]=1./9.; + w_lb[2]=1./9.; + w_lb[3]=1./9.; + w_lb[4]=1./9.; + w_lb[5]=1./9.; + w_lb[6]=1./9.; + w_lb[7]=1./72.; + w_lb[8]=1./72.; + w_lb[9]=1./72.; + w_lb[10]=1./72.; + w_lb[11]=1./72.; + w_lb[12]=1./72.; + w_lb[13]=1./72.; + w_lb[14]=1./72.; + + Ng_lb[0]=1.; + Ng_lb[1]=3.; + Ng_lb[2]=3.; + Ng_lb[3]=3.; + Ng_lb[4]=9./2.; + Ng_lb[5]=9./2.; + Ng_lb[6]=9./2.; + Ng_lb[7]=9.; + Ng_lb[8]=9.; + Ng_lb[9]=9.; + Ng_lb[10]=27./2.; + Ng_lb[11]=27./2.; + Ng_lb[12]=27./2.; + Ng_lb[13]=9.; + Ng_lb[14]=1.; + + mg_lb[0][0]=1.;mg_lb[0][1]=1.;mg_lb[0][2]=1.;mg_lb[0][3]=1.;mg_lb[0][4]=1.; + mg_lb[0][5]=1.;mg_lb[0][6]=1.;mg_lb[0][7]=1.;mg_lb[0][8]=1.;mg_lb[0][9]=1.; + mg_lb[0][10]=1.;mg_lb[0][11]=1.;mg_lb[0][12]=1.;mg_lb[0][13]=1.;mg_lb[0][14]=1.; + mg_lb[1][0]=0;mg_lb[1][1]=1.;mg_lb[1][2]=0;mg_lb[1][3]=-1.;mg_lb[1][4]=0; + mg_lb[1][5]=0;mg_lb[1][6]=0;mg_lb[1][7]=1.;mg_lb[1][8]=-1.;mg_lb[1][9]=-1.; + mg_lb[1][10]=1.;mg_lb[1][11]=1.;mg_lb[1][12]=-1.;mg_lb[1][13]=-1.;mg_lb[1][14]=1.; + mg_lb[2][0]=0;mg_lb[2][1]=0;mg_lb[2][2]=1.;mg_lb[2][3]=0;mg_lb[2][4]=-1.; + mg_lb[2][5]=0;mg_lb[2][6]=0;mg_lb[2][7]=1.;mg_lb[2][8]=1.;mg_lb[2][9]=-1.; + mg_lb[2][10]=-1.;mg_lb[2][11]=1.;mg_lb[2][12]=1.;mg_lb[2][13]=-1.;mg_lb[2][14]=-1.; + mg_lb[3][0]=0;mg_lb[3][1]=0;mg_lb[3][2]=0;mg_lb[3][3]=0;mg_lb[3][4]=0; + mg_lb[3][5]=1.;mg_lb[3][6]=-1.;mg_lb[3][7]=1.;mg_lb[3][8]=1.;mg_lb[3][9]=1.; + mg_lb[3][10]=1.;mg_lb[3][11]=-1.;mg_lb[3][12]=-1.;mg_lb[3][13]=-1.;mg_lb[3][14]=-1.; + mg_lb[4][0]=-1./3.;mg_lb[4][1]=2./3.;mg_lb[4][2]=-1./3.;mg_lb[4][3]=2./3.;mg_lb[4][4]=-1./3.; + mg_lb[4][5]=-1./3.;mg_lb[4][6]=-1./3.;mg_lb[4][7]=2./3.;mg_lb[4][8]=2./3.;mg_lb[4][9]=2./3.; + mg_lb[4][10]=2./3.;mg_lb[4][11]=2./3.;mg_lb[4][12]=2./3.;mg_lb[4][13]=2./3.;mg_lb[4][14]=2./3.; + mg_lb[5][0]=-1./3.;mg_lb[5][1]=-1./3.;mg_lb[5][2]=2./3.;mg_lb[5][3]=-1./3.;mg_lb[5][4]=2./3.; + mg_lb[5][5]=-1./3.;mg_lb[5][6]=-1./3.;mg_lb[5][7]=2./3.;mg_lb[5][8]=2./3.;mg_lb[5][9]=2./3.; + mg_lb[5][10]=2./3.;mg_lb[5][11]=2./3.;mg_lb[5][12]=2./3.;mg_lb[5][13]=2./3.;mg_lb[5][14]=2./3.; + mg_lb[6][0]=-1./3.;mg_lb[6][1]=-1./3.;mg_lb[6][2]=-1./3.;mg_lb[6][3]=-1./3.;mg_lb[6][4]=-1./3.; + mg_lb[6][5]=2./3.;mg_lb[6][6]=2./3.;mg_lb[6][7]=2./3.;mg_lb[6][8]=2./3.;mg_lb[6][9]=2./3.; + mg_lb[6][10]=2./3.;mg_lb[6][11]=2./3.;mg_lb[6][12]=2./3.;mg_lb[6][13]=2./3.;mg_lb[6][14]=2./3.; + mg_lb[7][0]=0;mg_lb[7][1]=0;mg_lb[7][2]=0;mg_lb[7][3]=0;mg_lb[7][4]=0; + mg_lb[7][5]=0;mg_lb[7][6]=0;mg_lb[7][7]=1;mg_lb[7][8]=-1;mg_lb[7][9]=1; + mg_lb[7][10]=-1;mg_lb[7][11]=1;mg_lb[7][12]=-1;mg_lb[7][13]=1;mg_lb[7][14]=-1; + mg_lb[8][0]=0;mg_lb[8][1]=0;mg_lb[8][2]=0;mg_lb[8][3]=0;mg_lb[8][4]=0; + mg_lb[8][5]=0;mg_lb[8][6]=0;mg_lb[8][7]=1;mg_lb[8][8]=1;mg_lb[8][9]=-1; + mg_lb[8][10]=-1;mg_lb[8][11]=-1;mg_lb[8][12]=-1;mg_lb[8][13]=1;mg_lb[8][14]=1; + mg_lb[9][0]=0;mg_lb[9][1]=0;mg_lb[9][2]=0;mg_lb[9][3]=0;mg_lb[9][4]=0; + mg_lb[9][5]=0;mg_lb[9][6]=0;mg_lb[9][7]=1;mg_lb[9][8]=-1;mg_lb[9][9]=-1; + mg_lb[9][10]=1;mg_lb[9][11]=-1;mg_lb[9][12]=1;mg_lb[9][13]=1;mg_lb[9][14]=-1; + mg_lb[10][0]=0;mg_lb[10][1]=0;mg_lb[10][2]=-1./3.;mg_lb[10][3]=0;mg_lb[10][4]=1./3.; + mg_lb[10][5]=0;mg_lb[10][6]=0;mg_lb[10][7]=2./3.;mg_lb[10][8]=2./3.;mg_lb[10][9]=-2./3.; + mg_lb[10][10]=-2./3.;mg_lb[10][11]=2./3.;mg_lb[10][12]=2./3.;mg_lb[10][13]=-2./3.;mg_lb[10][14]=-2./3.; + mg_lb[11][0]=0;mg_lb[11][1]=0;mg_lb[11][2]=0;mg_lb[11][3]=0;mg_lb[11][4]=0; + mg_lb[11][5]=-1./3.;mg_lb[11][6]=1./3.;mg_lb[11][7]=2./3.;mg_lb[11][8]=2./3.;mg_lb[11][9]=2./3.; + mg_lb[11][10]=2./3.;mg_lb[11][11]=-2./3.;mg_lb[11][12]=-2./3.;mg_lb[11][13]=-2./3.;mg_lb[11][14]=-2./3.; + mg_lb[12][0]=0;mg_lb[12][1]=-1./3.;mg_lb[12][2]=0;mg_lb[12][3]=1./3.;mg_lb[12][4]=0; + mg_lb[12][5]=0;mg_lb[12][6]=0;mg_lb[12][7]=2./3.;mg_lb[12][8]=-2./3.;mg_lb[12][9]=-2./3.; + mg_lb[12][10]=2./3.;mg_lb[12][11]=2./3.;mg_lb[12][12]=-2./3.;mg_lb[12][13]=-2./3.;mg_lb[12][14]=2./3.; + mg_lb[13][0]=0;mg_lb[13][1]=0;mg_lb[13][2]=0;mg_lb[13][3]=0;mg_lb[13][4]=0; + mg_lb[13][5]=0;mg_lb[13][6]=0;mg_lb[13][7]=1;mg_lb[13][8]=-1;mg_lb[13][9]=1; + mg_lb[13][10]=-1;mg_lb[13][11]=-1;mg_lb[13][12]=1;mg_lb[13][13]=-1;mg_lb[13][14]=1; + mg_lb[14][0]=sqrt(2.);mg_lb[14][1]=-1./sqrt(2.);mg_lb[14][2]=-1./sqrt(2.); + mg_lb[14][3]=-1./sqrt(2.);mg_lb[14][4]=-1./sqrt(2.); + mg_lb[14][5]=-1./sqrt(2.);mg_lb[14][6]=-1./sqrt(2.);mg_lb[14][7]=sqrt(2.); + mg_lb[14][8]=sqrt(2.);mg_lb[14][9]=sqrt(2.); + mg_lb[14][10]=sqrt(2.);mg_lb[14][11]=sqrt(2.);mg_lb[14][12]=sqrt(2.); + mg_lb[14][13]=sqrt(2.);mg_lb[14][14]=sqrt(2.); + + for(i=0; iprocneigh[0][0],15,world,&requests[0]); + MPI_Irecv(&feq[0][1][1][0],1,passxf,comm->procneigh[0][0],25,world,&requests[1]); + MPI_Isend(&feq[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],25,world,&requests[2]); + MPI_Irecv(&feq[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],15,world,&requests[3]); + if(typeLB == 2){ + MPI_Isend(&feqn[1][1][1][0],1,passxf,comm->procneigh[0][0],10,world,&requests[4]); + MPI_Irecv(&feqn[0][1][1][0],1,passxf,comm->procneigh[0][0],20,world,&requests[5]); + MPI_Isend(&feqn[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],20,world,&requests[6]); + MPI_Irecv(&feqn[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],10,world,&requests[7]); + } + MPI_Waitall(numrequests,requests,statuses); + + for(i=0; iprocneigh[1][0],15,world,&requests[0]); + MPI_Irecv(&feq[0][0][1][0],1,passyf,comm->procneigh[1][0],25,world,&requests[1]); + MPI_Isend(&feq[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],25,world,&requests[2]); + MPI_Irecv(&feq[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],15,world,&requests[3]); + if(typeLB == 2){ + MPI_Isend(&feqn[0][1][1][0],1,passyf,comm->procneigh[1][0],10,world,&requests[4]); + MPI_Irecv(&feqn[0][0][1][0],1,passyf,comm->procneigh[1][0],20,world,&requests[5]); + MPI_Isend(&feqn[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],20,world,&requests[6]); + MPI_Irecv(&feqn[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],10,world,&requests[7]); + } + MPI_Waitall(numrequests,requests,statuses); + + for(i=0; iprocneigh[2][0],15,world,&requests[0]); + MPI_Irecv(&feq[0][0][0][0],1,passzf,comm->procneigh[2][0],25,world,&requests[1]); + MPI_Isend(&feq[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],25,world,&requests[2]); + MPI_Irecv(&feq[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],15,world,&requests[3]); + if(typeLB == 2){ + MPI_Isend(&feqn[0][0][1][0],1,passzf,comm->procneigh[2][0],10,world,&requests[4]); + MPI_Irecv(&feqn[0][0][0][0],1,passzf,comm->procneigh[2][0],20,world,&requests[5]); + MPI_Isend(&feqn[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],20,world,&requests[6]); + MPI_Irecv(&feqn[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],10,world,&requests[7]); + } + MPI_Waitall(numrequests,requests,statuses); + + //Save feqold. + if(typeLB == 2){ + for(i=0; iprocneigh[0][0],15,world,&requests[0]); + MPI_Irecv(&feqold[0][1][1][0],1,passxf,comm->procneigh[0][0],25,world,&requests[1]); + MPI_Isend(&feqold[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],25,world,&requests[2]); + MPI_Irecv(&feqold[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],15,world,&requests[3]); + MPI_Isend(&feqoldn[1][1][1][0],1,passxf,comm->procneigh[0][0],10,world,&requests[4]); + MPI_Irecv(&feqoldn[0][1][1][0],1,passxf,comm->procneigh[0][0],20,world,&requests[5]); + MPI_Isend(&feqoldn[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],20,world,&requests[6]); + MPI_Irecv(&feqoldn[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],10,world,&requests[7]); + MPI_Waitall(8,requests,statuses); + + for(i=0; i<8; i++) + requests[i]=MPI_REQUEST_NULL; + MPI_Isend(&feqold[0][1][1][0],1,passyf,comm->procneigh[1][0],15,world,&requests[0]); + MPI_Irecv(&feqold[0][0][1][0],1,passyf,comm->procneigh[1][0],25,world,&requests[1]); + MPI_Isend(&feqold[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],25,world,&requests[2]); + MPI_Irecv(&feqold[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],15,world,&requests[3]); + MPI_Isend(&feqoldn[0][1][1][0],1,passyf,comm->procneigh[1][0],10,world,&requests[4]); + MPI_Irecv(&feqoldn[0][0][1][0],1,passyf,comm->procneigh[1][0],20,world,&requests[5]); + MPI_Isend(&feqoldn[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],20,world,&requests[6]); + MPI_Irecv(&feqoldn[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],10,world,&requests[7]); + MPI_Waitall(8,requests,statuses); + + for(i=0; i<8; i++) + requests[i]=MPI_REQUEST_NULL; + MPI_Isend(&feqold[0][0][1][0],1,passzf,comm->procneigh[2][0],15,world,&requests[0]); + MPI_Irecv(&feqold[0][0][0][0],1,passzf,comm->procneigh[2][0],25,world,&requests[1]); + MPI_Isend(&feqold[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],25,world,&requests[2]); + MPI_Irecv(&feqold[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],15,world,&requests[3]); + MPI_Isend(&feqoldn[0][0][1][0],1,passzf,comm->procneigh[2][0],10,world,&requests[4]); + MPI_Irecv(&feqoldn[0][0][0][0],1,passzf,comm->procneigh[2][0],20,world,&requests[5]); + MPI_Isend(&feqoldn[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],20,world,&requests[6]); + MPI_Irecv(&feqoldn[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],10,world,&requests[7]); + MPI_Waitall(8,requests,statuses); + } + parametercalc_full(); + } +} + +//========================================================================== +// Compute the lattice Boltzmann equilibrium distribution functions for +// the D3Q15 model. +//========================================================================== +void FixLbFluid::equilibriumdist15(int xstart, int xend, int ystart, int yend, int zstart, int zend) { + + double rho; + double usq; + int i, j, k, l, iup, idwn, jup, jdwn, kup, kdwn; + double Fx_w, Fy_w, Fz_w; + + double total_density(0.0); + double drhox, drhoy, drhoz, drhoxx, drhoyy, drhozz; + double Pxx, Pyy, Pzz, Pxy, Pxz, Pyz; + double grs, p0,TrP; + double dPdrho; + + double S[2][3],std; + int jj; + + double etacov[15],ghostnoise; + + + for (i=xstart; iperiodicity[2]==0){ + if(comm->myloc[2]==0 && k==1){ + drhoz = (-3.0*density_lb[i][j][k] + 4.0*density_lb[i][j][k+1] - + density_lb[i][j][k+2])/2.0; + drhozz = (-density_lb[i][j][k+3] + 4.0*density_lb[i][j][k+2] - + 5.0*density_lb[i][j][k+1] + 2.0*rho); + } + if(comm->myloc[2]==comm->procgrid[2]-1 && k==subNbz-2){ + drhoz = -(-3.0*density_lb[i][j][k] + 4.0*density_lb[i][j][k-1] - + density_lb[i][j][k-2])/2.0; + drhozz = (-density_lb[i][j][k-3] + 4.0*density_lb[i][j][k-2] - + 5.0*density_lb[i][j][k-1] + 2.0*rho); + } + } + + grs = drhox*drhox + drhoy*drhoy + drhoz*drhoz; + + p0 = rho*a_0-kappa_lb*rho*(drhoxx + drhoyy + drhozz); +// kappa_lb is the square gradient coeff in the pressure tensor + + dPdrho = a_0; //assuming here that kappa_lb = 0. + + + if(typeLB==1){ + Pxx = p0 + kappa_lb*(drhox*drhox - 0.5*grs)+(tau-0.5)*(1.0/3.0-dPdrho)* + (3.0*u_lb[i][j][k][0]*drhox+u_lb[i][j][k][1]*drhoy+u_lb[i][j][k][2]*drhoz); + Pyy = p0 + kappa_lb*(drhoy*drhoy - 0.5*grs)+(tau-0.5)*(1.0/3.0-dPdrho)* + (u_lb[i][j][k][0]*drhox+3.0*u_lb[i][j][k][1]*drhoy+u_lb[i][j][k][2]*drhoz); + Pzz = p0 + kappa_lb*(drhoz*drhoz - 0.5*grs)+(tau-0.5)*(1.0/3.0-dPdrho)* + (u_lb[i][j][k][0]*drhox+u_lb[i][j][k][1]*drhoy+3.0*u_lb[i][j][k][2]*drhoz); + Pxy = kappa_lb*drhox*drhoy+(tau-0.5)*(1.0/3.0-dPdrho)* + (u_lb[i][j][k][0]*drhoy+u_lb[i][j][k][1]*drhox); + Pxz = kappa_lb*drhox*drhoz+(tau-0.5)*(1.0/3.0-dPdrho)* + (u_lb[i][j][k][0]*drhoz+u_lb[i][j][k][2]*drhox); + Pyz = kappa_lb*drhoy*drhoz+(tau-0.5)*(1.0/3.0-dPdrho)* + (u_lb[i][j][k][1]*drhoz+u_lb[i][j][k][2]*drhoy); + }else if(typeLB==2){ + Pxx = p0 + kappa_lb*(drhox*drhox - 0.5*grs)+tau*(1.0/3.0-dPdrho)* + (3.0*u_lb[i][j][k][0]*drhox+u_lb[i][j][k][1]*drhoy+u_lb[i][j][k][2]*drhoz); + Pyy = p0 + kappa_lb*(drhoy*drhoy - 0.5*grs)+tau*(1.0/3.0-dPdrho)* + (u_lb[i][j][k][0]*drhox+3.0*u_lb[i][j][k][1]*drhoy+u_lb[i][j][k][2]*drhoz); + Pzz = p0 + kappa_lb*(drhoz*drhoz - 0.5*grs)+tau*(1.0/3.0-dPdrho)* + (u_lb[i][j][k][0]*drhox+u_lb[i][j][k][1]*drhoy+3.0*u_lb[i][j][k][2]*drhoz); + Pxy = kappa_lb*drhox*drhoy+tau*(1.0/3.0-dPdrho)* + (u_lb[i][j][k][0]*drhoy+u_lb[i][j][k][1]*drhox); + Pxz = kappa_lb*drhox*drhoz+tau*(1.0/3.0-dPdrho)* + (u_lb[i][j][k][0]*drhoz+u_lb[i][j][k][2]*drhox); + Pyz = kappa_lb*drhoy*drhoz+tau*(1.0/3.0-dPdrho)* + (u_lb[i][j][k][1]*drhoz+u_lb[i][j][k][2]*drhoy); + } + + Fx_w = Ff[i][j][k][0]; + Fy_w = Ff[i][j][k][1]; + Fz_w = Ff[i][j][k][2]; + + // this is Tr(P) + TrP = Pxx+Pyy+Pzz; + usq=u_lb[i][j][k][0]*u_lb[i][j][k][0]+u_lb[i][j][k][1]*u_lb[i][j][k][1]+ + u_lb[i][j][k][2]*u_lb[i][j][k][2]; + + etacov[0] = rho; + etacov[1] = rho*u_lb[i][j][k][0] + Fx_w*tau + rho*bodyforcex*tau; + etacov[2] = rho*u_lb[i][j][k][1] + Fy_w*tau + rho*bodyforcey*tau; + etacov[3] = rho*u_lb[i][j][k][2] + Fz_w*tau + rho*bodyforcez*tau; + + etacov[4] = Pxx + rho*u_lb[i][j][k][0]*u_lb[i][j][k][0] -rho/3. + + tau*(2.0*u_lb[i][j][k][0]*(Fx_w+rho*bodyforcex)); + etacov[5] = Pyy + rho*u_lb[i][j][k][1]*u_lb[i][j][k][1] -rho/3. + + tau*(2.0*u_lb[i][j][k][1]*(Fy_w+rho*bodyforcey)); + etacov[6] = Pzz + rho*u_lb[i][j][k][2]*u_lb[i][j][k][2] -rho/3. + + tau*(2.0*u_lb[i][j][k][2]*(Fz_w+rho*bodyforcez)); + etacov[7] = Pxy + rho*u_lb[i][j][k][0]*u_lb[i][j][k][1] + + tau*(u_lb[i][j][k][0]*(Fy_w+rho*bodyforcey) + (Fx_w+rho*bodyforcex)*u_lb[i][j][k][1]); + etacov[8] = Pyz + rho*u_lb[i][j][k][1]*u_lb[i][j][k][2] + + tau*(u_lb[i][j][k][1]*(Fz_w+rho*bodyforcez) + (Fy_w+rho*bodyforcey)*u_lb[i][j][k][2]); + etacov[9] = Pxz + rho*u_lb[i][j][k][0]*u_lb[i][j][k][2] + + tau*(u_lb[i][j][k][0]*(Fz_w+rho*bodyforcez) + (Fx_w+rho*bodyforcex)*u_lb[i][j][k][2]); + 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]*u_lb[i][j][k][2]; + etacov[14] = K_0*(rho-TrP); + + for (l=0; l<15; l++) { + + feq[i][j][k][l] = 0.0; + for (int ii=0; ii<15; ii++) + feq[i][j][k][l] += w_lb[l]*mg_lb[ii][l]*etacov[ii]*Ng_lb[ii]; + + if(typeLB == 2){ + feqn[i][j][k][l] = feq[i][j][k][l]; + } + } + + if(noisestress==1){ + std = sqrt(namp*rho); + + for(jj=0; jj<3; jj++) + S[0][jj] = std*random->gaussian(); + for(jj=0; jj<3; jj++) + S[1][jj] = std*random->gaussian(); + + etacov[4] = (S[0][0]*sqrt(3.0-3.0*a_0)); + etacov[5] = ((1.0-3.0*a_0)*S[0][0]/sqrt(3.0-3.0*a_0)+ + sqrt((8.0-12.0*a_0)/(3.0-3.0*a_0))*S[0][1]); + etacov[6] = ((1.0-3.0*a_0)*S[0][0]/sqrt(3.0-3.0*a_0)+ + (2.0-6.0*a_0)*S[0][1]/sqrt((8.0-12.0*a_0)*(3.0-3.0*a_0))+ + sqrt((5.0-9.0*a_0)/(2.0-3.0*a_0))*S[0][2]); + etacov[7] = S[1][0]; + etacov[8] = S[1][1]; + etacov[9] = S[1][2]; + + for (l=10; l<15; l++) { + etacov[l] = sqrt(9.0*namp*rho/Ng_lb[l])*random->gaussian(); + } + etacov[14] += -K_0*(etacov[4]+etacov[5]+etacov[6]); //correction from noise to TrP + + for (l=0; l<15; l++) { + ghostnoise = w_lb[l]* + (mg_lb[4][l]*etacov[4]*Ng_lb[4] + mg_lb[5][l]*etacov[5]*Ng_lb[5] + + mg_lb[6][l]*etacov[6]*Ng_lb[6] + mg_lb[7][l]*etacov[7]*Ng_lb[7] + + mg_lb[8][l]*etacov[8]*Ng_lb[8] + mg_lb[9][l]*etacov[9]*Ng_lb[9] + + mg_lb[10][l]*etacov[10]*Ng_lb[10] + mg_lb[11][l]*etacov[11]*Ng_lb[11] + + mg_lb[12][l]*etacov[12]*Ng_lb[12] + mg_lb[13][l]*etacov[13]*Ng_lb[13] + + mg_lb[14][l]*etacov[14]*Ng_lb[14]); + feq[i][j][k][l] += ghostnoise*noisefactor; + } + } + } + } + } +} + +//========================================================================== +// Compute the lattice Boltzmann equilibrium distribution functions for +// the D3Q19 model. +//========================================================================== +void FixLbFluid::equilibriumdist19(int xstart, int xend, int ystart, int yend, int zstart, int zend) { + + double rho; + double usq; + int i, j, k, l, iup, idwn, jup, jdwn, kup, kdwn; + double Fx_w, Fy_w, Fz_w; + + double total_density(0.0); + double drhox, drhoy, drhoz, drhoxx, drhoyy, drhozz; + double Pxx, Pyy, Pzz, Pxy, Pxz, Pyz; + double grs, p0,TrP; + double dPdrho; + + double S[2][3],std; + int jj; + + double etacov[19],ghostnoise; + + for (i=xstart; iperiodicity[2]==0){ + if(comm->myloc[2]==0 && k==1){ + drhoz = (-3.0*density_lb[i][j][k] + 4.0*density_lb[i][j][k+1] - + density_lb[i][j][k+2])/2.0; + drhozz = (-density_lb[i][j][k+3] + 4.0*density_lb[i][j][k+2] - + 5.0*density_lb[i][j][k+1] + 2.0*rho); + } + if(comm->myloc[2]==comm->procgrid[2]-1 && k==subNbz-2){ + drhoz = -(-3.0*density_lb[i][j][k] + 4.0*density_lb[i][j][k-1] - + density_lb[i][j][k-2])/2.0; + drhozz = (-density_lb[i][j][k-3] + 4.0*density_lb[i][j][k-2] - + 5.0*density_lb[i][j][k-1] + 2.0*rho); + } + } + + grs = drhox*drhox + drhoy*drhoy + drhoz*drhoz; + + p0 = rho*a_0-kappa_lb*rho*(drhoxx + drhoyy + drhozz); +// kappa_lb is the square gradient coeff in the pressure tensor + + dPdrho = a_0; //assuming here that kappa_lb = 0. + + + if(typeLB==1){ + Pxx = p0 + kappa_lb*(drhox*drhox - 0.5*grs)+(tau-0.5)*(1.0/3.0-dPdrho)* + (3.0*u_lb[i][j][k][0]*drhox+u_lb[i][j][k][1]*drhoy+u_lb[i][j][k][2]*drhoz); + Pyy = p0 + kappa_lb*(drhoy*drhoy - 0.5*grs)+(tau-0.5)*(1.0/3.0-dPdrho)* + (u_lb[i][j][k][0]*drhox+3.0*u_lb[i][j][k][1]*drhoy+u_lb[i][j][k][2]*drhoz); + Pzz = p0 + kappa_lb*(drhoz*drhoz - 0.5*grs)+(tau-0.5)*(1.0/3.0-dPdrho)* + (u_lb[i][j][k][0]*drhox+u_lb[i][j][k][1]*drhoy+3.0*u_lb[i][j][k][2]*drhoz); + Pxy = kappa_lb*drhox*drhoy+(tau-0.5)*(1.0/3.0-dPdrho)* + (u_lb[i][j][k][0]*drhoy+u_lb[i][j][k][1]*drhox); + Pxz = kappa_lb*drhox*drhoz+(tau-0.5)*(1.0/3.0-dPdrho)* + (u_lb[i][j][k][0]*drhoz+u_lb[i][j][k][2]*drhox); + Pyz = kappa_lb*drhoy*drhoz+(tau-0.5)*(1.0/3.0-dPdrho)* + (u_lb[i][j][k][1]*drhoz+u_lb[i][j][k][2]*drhoy); + }else if(typeLB==2){ + Pxx = p0 + kappa_lb*(drhox*drhox - 0.5*grs)+tau*(1.0/3.0-dPdrho)* + (3.0*u_lb[i][j][k][0]*drhox+u_lb[i][j][k][1]*drhoy+u_lb[i][j][k][2]*drhoz); + Pyy = p0 + kappa_lb*(drhoy*drhoy - 0.5*grs)+tau*(1.0/3.0-dPdrho)* + (u_lb[i][j][k][0]*drhox+3.0*u_lb[i][j][k][1]*drhoy+u_lb[i][j][k][2]*drhoz); + Pzz = p0 + kappa_lb*(drhoz*drhoz - 0.5*grs)+tau*(1.0/3.0-dPdrho)* + (u_lb[i][j][k][0]*drhox+u_lb[i][j][k][1]*drhoy+3.0*u_lb[i][j][k][2]*drhoz); + Pxy = kappa_lb*drhox*drhoy+tau*(1.0/3.0-dPdrho)* + (u_lb[i][j][k][0]*drhoy+u_lb[i][j][k][1]*drhox); + Pxz = kappa_lb*drhox*drhoz+tau*(1.0/3.0-dPdrho)* + (u_lb[i][j][k][0]*drhoz+u_lb[i][j][k][2]*drhox); + Pyz = kappa_lb*drhoy*drhoz+tau*(1.0/3.0-dPdrho)* + (u_lb[i][j][k][1]*drhoz+u_lb[i][j][k][2]*drhoy); + } + + Fx_w = Ff[i][j][k][0]; + Fy_w = Ff[i][j][k][1]; + Fz_w = Ff[i][j][k][2]; + + // this is Tr(P) + TrP = Pxx+Pyy+Pzz; + usq=u_lb[i][j][k][0]*u_lb[i][j][k][0]+u_lb[i][j][k][1]*u_lb[i][j][k][1]+ + u_lb[i][j][k][2]*u_lb[i][j][k][2]; + + etacov[0] = rho; + etacov[1] = rho*u_lb[i][j][k][0] + Fx_w*tau + rho*bodyforcex*tau; + etacov[2] = rho*u_lb[i][j][k][1] + Fy_w*tau + rho*bodyforcey*tau; + etacov[3] = rho*u_lb[i][j][k][2] + Fz_w*tau + rho*bodyforcez*tau; + + etacov[4] = Pxx + rho*u_lb[i][j][k][0]*u_lb[i][j][k][0] -rho/3. + + tau*(2.0*u_lb[i][j][k][0]*(Fx_w+rho*bodyforcex)); + etacov[5] = Pyy + rho*u_lb[i][j][k][1]*u_lb[i][j][k][1] -rho/3. + + tau*(2.0*u_lb[i][j][k][1]*(Fy_w+rho*bodyforcey)); + etacov[6] = Pzz + rho*u_lb[i][j][k][2]*u_lb[i][j][k][2] -rho/3. + + tau*(2.0*u_lb[i][j][k][2]*(Fz_w+rho*bodyforcez)); + etacov[7] = Pxy + rho*u_lb[i][j][k][0]*u_lb[i][j][k][1] + + tau*(u_lb[i][j][k][0]*(Fy_w+rho*bodyforcey) + (Fx_w+rho*bodyforcex)*u_lb[i][j][k][1]); + etacov[8] = Pxz + rho*u_lb[i][j][k][0]*u_lb[i][j][k][2] + + tau*(u_lb[i][j][k][0]*(Fz_w+rho*bodyforcez) + (Fx_w+rho*bodyforcex)*u_lb[i][j][k][2]); + etacov[9] = Pyz + rho*u_lb[i][j][k][1]*u_lb[i][j][k][2] + + tau*(u_lb[i][j][k][1]*(Fz_w+rho*bodyforcez) + (Fy_w+rho*bodyforcey)*u_lb[i][j][k][2]); + 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 (l=0; l<19; l++) { + + feq[i][j][k][l] = 0.0; + for (int ii=0; ii<19; ii++) + feq[i][j][k][l] += w_lb[l]*mg_lb[ii][l]*etacov[ii]*Ng_lb[ii]; + + if(typeLB == 2){ + feqn[i][j][k][l] = feq[i][j][k][l]; + } + } + + if(noisestress==1){ + std = sqrt(namp*rho); + + for(jj=0; jj<3; jj++) + S[0][jj] = std*random->gaussian(); + for(jj=0; jj<3; jj++) + S[1][jj] = std*random->gaussian(); + + etacov[4] = (S[0][0]*sqrt(3.0-3.0*a_0)); + etacov[5] = ((1.0-3.0*a_0)*S[0][0]/sqrt(3.0-3.0*a_0)+ + sqrt((8.0-12.0*a_0)/(3.0-3.0*a_0))*S[0][1]); + etacov[6] = ((1.0-3.0*a_0)*S[0][0]/sqrt(3.0-3.0*a_0)+ + (2.0-6.0*a_0)*S[0][1]/sqrt((8.0-12.0*a_0)*(3.0-3.0*a_0))+ + sqrt((5.0-9.0*a_0)/(2.0-3.0*a_0))*S[0][2]); + etacov[7] = S[1][0]; + etacov[8] = S[1][1]; + etacov[9] = S[1][2]; + + for (l=10; l<19; l++) { + etacov[l] = sqrt(9.0*namp*rho/Ng_lb[l])*random->gaussian(); + } + + for (l=0; l<19; l++) { + ghostnoise = w_lb[l]* + (mg_lb[4][l]*etacov[4]*Ng_lb[4] + mg_lb[5][l]*etacov[5]*Ng_lb[5] + + mg_lb[6][l]*etacov[6]*Ng_lb[6] + mg_lb[7][l]*etacov[7]*Ng_lb[7] + + mg_lb[8][l]*etacov[8]*Ng_lb[8] + mg_lb[9][l]*etacov[9]*Ng_lb[9] + + mg_lb[10][l]*etacov[10]*Ng_lb[10] + mg_lb[11][l]*etacov[11]*Ng_lb[11] + + mg_lb[12][l]*etacov[12]*Ng_lb[12] + mg_lb[13][l]*etacov[13]*Ng_lb[13] + + mg_lb[14][l]*etacov[14]*Ng_lb[14] + mg_lb[15][l]*etacov[15]*Ng_lb[15] + + mg_lb[16][l]*etacov[16]*Ng_lb[16] + mg_lb[17][l]*etacov[17]*Ng_lb[17] + + mg_lb[18][l]*etacov[18]*Ng_lb[18]); + feq[i][j][k][l] += ghostnoise*noisefactor; + } + } + + } + } + } +} + +//========================================================================== +// Calculate the fluid density and velocity over the entire simulation +// domain. +//========================================================================== +void FixLbFluid::parametercalc_full(void) +{ + MPI_Request requests[4]; + MPI_Status statuses[4]; + MPI_Request requests2[12]; + MPI_Status statuses2[12]; + int numrequests; + int i; + + //-------------------------------------------------------------------------- + // send the boundaries of f_lb, as they will be needed later by the update + // routine, and use these to calculate the density and velocity on the + // boundary. + //-------------------------------------------------------------------------- + for(i=0; i<4; i++) + requests[i]=MPI_REQUEST_NULL; + MPI_Isend(&f_lb[1][1][1][0],1,passxf,comm->procneigh[0][0],10,world,&requests[0]); + MPI_Irecv(&f_lb[0][1][1][0],1,passxf,comm->procneigh[0][0],20,world,&requests[1]); + MPI_Isend(&f_lb[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],20,world,&requests[2]); + MPI_Irecv(&f_lb[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],10,world,&requests[3]); + parametercalc_part(1,subNbx-1,1,subNby-1,1,subNbz-1); + MPI_Waitall(4,requests,statuses); + + for(i=0; i<4; i++) + requests[i]=MPI_REQUEST_NULL; + MPI_Isend(&f_lb[0][1][1][0],1,passyf,comm->procneigh[1][0],10,world,&requests[0]); + MPI_Irecv(&f_lb[0][0][1][0],1,passyf,comm->procneigh[1][0],20,world,&requests[1]); + MPI_Isend(&f_lb[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],20,world,&requests[2]); + MPI_Irecv(&f_lb[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],10,world,&requests[3]); + parametercalc_part(0,1,1,subNby-1,1,subNbz-1); + parametercalc_part(subNbx-1,subNbx,1,subNby-1,1,subNbz-1); + MPI_Waitall(4,requests,statuses); + + for(i=0; i<4; i++) + requests[i]=MPI_REQUEST_NULL; + MPI_Isend(&f_lb[0][0][1][0],1,passzf,comm->procneigh[2][0],10,world,&requests[0]); + MPI_Irecv(&f_lb[0][0][0][0],1,passzf,comm->procneigh[2][0],20,world,&requests[1]); + MPI_Isend(&f_lb[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],20,world,&requests[2]); + MPI_Irecv(&f_lb[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],10,world,&requests[3]); + parametercalc_part(0,subNbx,0,1,1,subNbz-1); + parametercalc_part(0,subNbx,subNby-1,subNby,1,subNbz-1); + MPI_Waitall(4,requests,statuses); + + parametercalc_part(0,subNbx,0,subNby,0,1); + parametercalc_part(0,subNbx,0,subNby,subNbz-1,subNbz); + + //-------------------------------------------------------------------------- + // Send the remaining portions of the u array (and density array if Gamma + // is set the default way). + //-------------------------------------------------------------------------- + if(setGamma == 0) numrequests = 12; + else numrequests = 6; + + for(i=0; iprocneigh[0][0],10,world,&requests2[0]); + MPI_Isend(&u_lb[3][0][0][0],1,passxu,comm->procneigh[0][0],20,world,&requests2[1]); + MPI_Isend(&u_lb[subNbx-3][0][0][0],1,passxu,comm->procneigh[0][1],30,world,&requests2[2]); + MPI_Irecv(&u_lb[subNbx][0][0][0],1,passxu,comm->procneigh[0][1],10,world,&requests2[3]); + MPI_Irecv(&u_lb[subNbx+1][0][0][0],1,passxu,comm->procneigh[0][1],20,world,&requests2[4]); + MPI_Irecv(&u_lb[subNbx+2][0][0][0],1,passxu,comm->procneigh[0][0],30,world,&requests2[5]); + if(setGamma==0){ + MPI_Isend(&density_lb[2][0][0],1,passxrho,comm->procneigh[0][0],40,world,&requests2[6]); + MPI_Isend(&density_lb[3][0][0],1,passxrho,comm->procneigh[0][0],50,world,&requests2[7]); + MPI_Isend(&density_lb[subNbx-3][0][0],1,passxrho,comm->procneigh[0][1],60,world,&requests2[8]); + MPI_Irecv(&density_lb[subNbx][0][0],1,passxrho,comm->procneigh[0][1],40,world,&requests2[9]); + MPI_Irecv(&density_lb[subNbx+1][0][0],1,passxrho,comm->procneigh[0][1],50,world,&requests2[10]); + MPI_Irecv(&density_lb[subNbx+2][0][0],1,passxrho,comm->procneigh[0][0],60,world,&requests2[11]); + } + MPI_Waitall(numrequests,requests2,statuses2); + + for(i=0; iprocneigh[1][0],10,world,&requests2[0]); + MPI_Isend(&u_lb[0][3][0][0],1,passyu,comm->procneigh[1][0],20,world,&requests2[1]); + MPI_Isend(&u_lb[0][subNby-3][0][0],1,passyu,comm->procneigh[1][1],30,world,&requests2[2]); + MPI_Irecv(&u_lb[0][subNby][0][0],1,passyu,comm->procneigh[1][1],10,world,&requests2[3]); + MPI_Irecv(&u_lb[0][subNby+1][0][0],1,passyu,comm->procneigh[1][1],20,world,&requests2[4]); + MPI_Irecv(&u_lb[0][subNby+2][0][0],1,passyu,comm->procneigh[1][0],30,world,&requests2[5]); + if(setGamma==0){ + MPI_Isend(&density_lb[0][2][0],1,passyrho,comm->procneigh[1][0],40,world,&requests2[6]); + MPI_Isend(&density_lb[0][3][0],1,passyrho,comm->procneigh[1][0],50,world,&requests2[7]); + MPI_Isend(&density_lb[0][subNby-3][0],1,passyrho,comm->procneigh[1][1],60,world,&requests2[8]); + MPI_Irecv(&density_lb[0][subNby][0],1,passyrho,comm->procneigh[1][1],40,world,&requests2[9]); + MPI_Irecv(&density_lb[0][subNby+1][0],1,passyrho,comm->procneigh[1][1],50,world,&requests2[10]); + MPI_Irecv(&density_lb[0][subNby+2][0],1,passyrho,comm->procneigh[1][0],60,world,&requests2[11]); + } + MPI_Waitall(numrequests,requests2,statuses2); + + for(i=0; i<12; i++) + requests2[i]=MPI_REQUEST_NULL; + int requestcount=0; + if(domain->periodicity[2]!=0 || comm->myloc[2] != 0){ + MPI_Isend(&u_lb[0][0][2][0],1,passzu,comm->procneigh[2][0],10,world,&requests2[requestcount]); + MPI_Isend(&u_lb[0][0][3][0],1,passzu,comm->procneigh[2][0],20,world,&requests2[requestcount+1]); + MPI_Irecv(&u_lb[0][0][subNbz+2][0],1,passzu,comm->procneigh[2][0],30,world,&requests2[requestcount+2]); + requestcount=requestcount+3; + if(setGamma==0){ + MPI_Isend(&density_lb[0][0][2],1,passzrho,comm->procneigh[2][0],40,world,&requests2[requestcount]); + MPI_Isend(&density_lb[0][0][3],1,passzrho,comm->procneigh[2][0],50,world,&requests2[requestcount+1]); + MPI_Irecv(&density_lb[0][0][subNbz+2],1,passzrho,comm->procneigh[2][0],60,world,&requests2[requestcount+2]); + requestcount=requestcount+3; + } + } + if(domain->periodicity[2]!=0 || comm->myloc[2] != (comm->procgrid[2]-1)){ + MPI_Isend(&u_lb[0][0][subNbz-3][0],1,passzu,comm->procneigh[2][1],30,world,&requests2[requestcount]); + MPI_Irecv(&u_lb[0][0][subNbz][0],1,passzu,comm->procneigh[2][1],10,world,&requests2[requestcount+1]); + MPI_Irecv(&u_lb[0][0][subNbz+1][0],1,passzu,comm->procneigh[2][1],20,world,&requests2[requestcount+2]); + requestcount=requestcount+3; + if(setGamma==0){ + MPI_Isend(&density_lb[0][0][subNbz-3],1,passzrho,comm->procneigh[2][1],60,world,&requests2[requestcount]); + MPI_Irecv(&density_lb[0][0][subNbz],1,passzrho,comm->procneigh[2][1],40,world,&requests2[requestcount+1]); + MPI_Irecv(&density_lb[0][0][subNbz+1],1,passzrho,comm->procneigh[2][1],50,world,&requests2[requestcount+2]); + requestcount=requestcount+3; + } + } + MPI_Waitall(requestcount,requests2,statuses2); + +} + +//========================================================================== +// Calculate the fluid density and velocity over a simulation volume +// specified by xstart,xend; ystart,yend; zstart,zend. +//========================================================================== +void FixLbFluid::parametercalc_part(int xstart, int xend, int ystart, int yend, int zstart, int zend) +{ + int i,j,k,m; + + for(i=xstart; iperiodicity[2]==0){ + if(comm->myloc[2]==0){ + if(k==1){ + u_lb[i][j][k][2]=0.0; + } + } + if(comm->myloc[2]==comm->procgrid[2]-1){ + if(k==subNbz-2){ + u_lb[i][j][k][2]=0.0; + } + } + + } + + u_lb[i][j][k][0]=u_lb[i][j][k][0]/density_lb[i][j][k]; + u_lb[i][j][k][1]=u_lb[i][j][k][1]/density_lb[i][j][k]; + u_lb[i][j][k][2]=u_lb[i][j][k][2]/density_lb[i][j][k]; + } + } + } + +} + +//========================================================================== +// Update the distribution function over a simulation volume specified +// by xstart,xend; ystart,yend; zstart,zend. +//========================================================================== +void FixLbFluid::update_periodic(int xstart, int xend, int ystart, int yend, int zstart, int zend) +{ + int i,j,k,m; + int imod,jmod,kmod,imodm,jmodm,kmodm; + + for(i=xstart; ime==0){ + // printf("%16.12f %16.12f %16.12f %16.12f\n",mass*dm_lb,momentum[0]*dm_lb*dx_lb/dt_lb,momentum[1]*dm_lb*dx_lb/dt_lb,momentum[2]*dm_lb*dx_lb/dt_lb); + // } + + sizeloc=(subNbx*subNby*subNbz*4); + MPI_Allreduce(&sizeloc,&size,1,MPI_INT,MPI_MAX,world); + + if(me==0){ + for(iproc=0; iproc < comm->nprocs; iproc++){ + if(iproc){ + MPI_Irecv(&buf[0][0][0][0],size,MPI_DOUBLE,iproc,0,world,&request_recv); + MPI_Wait(&request_recv,&status); + + istart=static_cast (buf[0][0][0][0]); + jstart=static_cast (buf[0][0][0][1]); + kstart=static_cast (buf[0][0][0][2]); + iend=static_cast (buf[0][0][1][0]); + jend=static_cast (buf[0][0][1][1]); + kend=static_cast (buf[0][0][1][2]); + + for(i=istart; imyloc[0]*(subNbx-2); + jstart=comm->myloc[1]*(subNby-2); + if(domain->periodicity[2]==0){ + if(comm->myloc[2]==comm->procgrid[2]-1){ + kstart=comm->myloc[2]*(subNbz-3); + }else{ + kstart=comm->myloc[2]*(subNbz-2); + } + }else{ + kstart=comm->myloc[2]*(subNbz-2); + } + iend=istart+subNbx-2; + jend=jstart+subNby-2; + kend=kstart+subNbz-2; + for(i=0; iperiodicity[2]==0){ + + for(i=0; iprocneigh[0][0],15,world,&requests[0]); + MPI_Irecv(&feq[0][1][1][0],1,passxf,comm->procneigh[0][0],25,world,&requests[1]); + MPI_Isend(&feq[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],25,world,&requests[2]); + MPI_Irecv(&feq[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],15,world,&requests[3]); + if(typeLB == 2){ + MPI_Isend(&feqn[1][1][1][0],1,passxf,comm->procneigh[0][0],10,world,&requests[4]); + MPI_Irecv(&feqn[0][1][1][0],1,passxf,comm->procneigh[0][0],20,world,&requests[5]); + MPI_Isend(&feqn[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],20,world,&requests[6]); + MPI_Irecv(&feqn[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],10,world,&requests[7]); + } + update_periodic(2,subNbx-2,2,subNby-2,2,subNbz-2); + MPI_Waitall(numrequests,requests,statuses); + + + for(i=0; iprocneigh[1][0],15,world,&requests[0]); + MPI_Irecv(&feq[0][0][1][0],1,passyf,comm->procneigh[1][0],25,world,&requests[1]); + MPI_Isend(&feq[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],25,world,&requests[2]); + MPI_Irecv(&feq[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],15,world,&requests[3]); + if(typeLB == 2){ + MPI_Isend(&feqn[0][1][1][0],1,passyf,comm->procneigh[1][0],10,world,&requests[4]); + MPI_Irecv(&feqn[0][0][1][0],1,passyf,comm->procneigh[1][0],20,world,&requests[5]); + MPI_Isend(&feqn[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],20,world,&requests[6]); + MPI_Irecv(&feqn[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],10,world,&requests[7]); + } + update_periodic(1,2,2,subNby-2,2,subNbz-2); + update_periodic(subNbx-2,subNbx-1,2,subNby-2,2,subNbz-2); + MPI_Waitall(numrequests,requests,statuses); + + for(i=0; iprocneigh[2][0],15,world,&requests[0]); + MPI_Irecv(&feq[0][0][0][0],1,passzf,comm->procneigh[2][0],25,world,&requests[1]); + MPI_Isend(&feq[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],25,world,&requests[2]); + MPI_Irecv(&feq[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],15,world,&requests[3]); + if(typeLB == 2){ + MPI_Isend(&feqn[0][0][1][0],1,passzf,comm->procneigh[2][0],10,world,&requests[4]); + MPI_Irecv(&feqn[0][0][0][0],1,passzf,comm->procneigh[2][0],20,world,&requests[5]); + MPI_Isend(&feqn[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],20,world,&requests[6]); + MPI_Irecv(&feqn[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],10,world,&requests[7]); + } + update_periodic(1,subNbx-1,1,2,2,subNbz-2); + update_periodic(1,subNbx-1,subNby-2,subNby-1,2,subNbz-2); + MPI_Waitall(numrequests,requests,statuses); + + if(typeLB==1){ + update_periodic(1,subNbx-1,1,subNby-1,1,2); + update_periodic(1,subNbx-1,1,subNby-1,subNbz-2,subNbz-1); + }else if(typeLB==2){ + if(comm->myloc[2]==0){ + for(i=1; imyloc[2]==comm->procgrid[2]-1){ + for(i=1;imyloc[2]==0){ + MPI_Isend(&fnew[0][0][1][0],1,passzf,comm->procneigh[2][0],15,world,&req_send15); + MPI_Irecv(&fnew[0][0][0][0],1,passzf,comm->procneigh[2][0],25,world,&req_recv25); + } + + if(comm->myloc[2]==comm->procgrid[2]-1){ + MPI_Isend(&fnew[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],25,world,&req_send25); + MPI_Irecv(&fnew[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],15,world,&req_recv15); + } + if(comm->myloc[2]==0){ + MPI_Wait(&req_send15,&status); + MPI_Wait(&req_recv25,&status); + + for(i=1;imyloc[2]==comm->procgrid[2]-1){ + MPI_Wait(&req_send25,&status); + MPI_Wait(&req_recv15,&status); + + for(i=1;iprocneigh[0][0],15,world,&requests[0]); + MPI_Irecv(&feq[0][1][1][0],1,passxf,comm->procneigh[0][0],25,world,&requests[1]); + MPI_Isend(&feq[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],25,world,&requests[2]); + MPI_Irecv(&feq[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],15,world,&requests[3]); + if(typeLB == 2){ + MPI_Isend(&feqn[1][1][1][0],1,passxf,comm->procneigh[0][0],10,world,&requests[4]); + MPI_Irecv(&feqn[0][1][1][0],1,passxf,comm->procneigh[0][0],20,world,&requests[5]); + MPI_Isend(&feqn[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],20,world,&requests[6]); + MPI_Irecv(&feqn[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],10,world,&requests[7]); + } + update_periodic(2,subNbx-2,2,subNby-2,2,subNbz-2); + MPI_Waitall(numrequests,requests,statuses); + + for(i=0; iprocneigh[1][0],15,world,&requests[0]); + MPI_Irecv(&feq[0][0][1][0],1,passyf,comm->procneigh[1][0],25,world,&requests[1]); + MPI_Isend(&feq[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],25,world,&requests[2]); + MPI_Irecv(&feq[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],15,world,&requests[3]); + if(typeLB == 2){ + MPI_Isend(&feqn[0][1][1][0],1,passyf,comm->procneigh[1][0],10,world,&requests[4]); + MPI_Irecv(&feqn[0][0][1][0],1,passyf,comm->procneigh[1][0],20,world,&requests[5]); + MPI_Isend(&feqn[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],20,world,&requests[6]); + MPI_Irecv(&feqn[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],10,world,&requests[7]); + } + update_periodic(1,2,2,subNby-2,2,subNbz-2); + update_periodic(subNbx-2,subNbx-1,2,subNby-2,2,subNbz-2); + MPI_Waitall(numrequests,requests,statuses); + + for(i=0; iprocneigh[2][0],15,world,&requests[0]); + MPI_Irecv(&feq[0][0][0][0],1,passzf,comm->procneigh[2][0],25,world,&requests[1]); + MPI_Isend(&feq[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],25,world,&requests[2]); + MPI_Irecv(&feq[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],15,world,&requests[3]); + if(typeLB == 2){ + MPI_Isend(&feqn[0][0][1][0],1,passzf,comm->procneigh[2][0],10,world,&requests[4]); + MPI_Irecv(&feqn[0][0][0][0],1,passzf,comm->procneigh[2][0],20,world,&requests[5]); + MPI_Isend(&feqn[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],20,world,&requests[6]); + MPI_Irecv(&feqn[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],10,world,&requests[7]); + } + update_periodic(1,subNbx-1,1,2,2,subNbz-2); + update_periodic(1,subNbx-1,subNby-2,subNby-1,2,subNbz-2); + MPI_Waitall(numrequests,requests,statuses); + + update_periodic(1,subNbx-1,1,subNby-1,1,2); + update_periodic(1,subNbx-1,1,subNby-1,subNbz-2,subNbz-1); + } + +} + +//========================================================================== +// Update the distribution functions over the entire simulation domain for +// the D3Q19 model. +//========================================================================== +void FixLbFluid::update_full19(void) +{ + + MPI_Request req_send15,req_recv15; + MPI_Request req_send25,req_recv25; + MPI_Request requests[8]; + MPI_Status statuses[8]; + int numrequests; + double tmp1,tmp2,tmp3; + MPI_Status status; + double rb; + int i,j,k,m; + int imod,jmod,kmod; + int imodm,jmodm,kmodm; + + //-------------------------------------------------------------------------- + // If using the standard LB integrator, do not need to send info about feqn. + //-------------------------------------------------------------------------- + if(typeLB == 1){ + numrequests = 4; + }else{ + numrequests = 8; + } + + //-------------------------------------------------------------------------- + // Fixed z boundary conditions. + //-------------------------------------------------------------------------- + if(domain->periodicity[2]==0){ + + for(i=0; iprocneigh[0][0],15,world,&requests[0]); + MPI_Irecv(&feq[0][1][1][0],1,passxf,comm->procneigh[0][0],25,world,&requests[1]); + MPI_Isend(&feq[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],25,world,&requests[2]); + MPI_Irecv(&feq[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],15,world,&requests[3]); + if(typeLB == 2){ + MPI_Isend(&feqn[1][1][1][0],1,passxf,comm->procneigh[0][0],10,world,&requests[4]); + MPI_Irecv(&feqn[0][1][1][0],1,passxf,comm->procneigh[0][0],20,world,&requests[5]); + MPI_Isend(&feqn[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],20,world,&requests[6]); + MPI_Irecv(&feqn[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],10,world,&requests[7]); + } + update_periodic(2,subNbx-2,2,subNby-2,2,subNbz-2); + MPI_Waitall(numrequests,requests,statuses); + + for(i=0; iprocneigh[1][0],15,world,&requests[0]); + MPI_Irecv(&feq[0][0][1][0],1,passyf,comm->procneigh[1][0],25,world,&requests[1]); + MPI_Isend(&feq[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],25,world,&requests[2]); + MPI_Irecv(&feq[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],15,world,&requests[3]); + if(typeLB == 2){ + MPI_Isend(&feqn[0][1][1][0],1,passyf,comm->procneigh[1][0],10,world,&requests[4]); + MPI_Irecv(&feqn[0][0][1][0],1,passyf,comm->procneigh[1][0],20,world,&requests[5]); + MPI_Isend(&feqn[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],20,world,&requests[6]); + MPI_Irecv(&feqn[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],10,world,&requests[7]); + } + update_periodic(1,2,2,subNby-2,2,subNbz-2); + update_periodic(subNbx-2,subNbx-1,2,subNby-2,2,subNbz-2); + MPI_Waitall(numrequests,requests,statuses); + + for(i=0; iprocneigh[2][0],15,world,&requests[0]); + MPI_Irecv(&feq[0][0][0][0],1,passzf,comm->procneigh[2][0],25,world,&requests[1]); + MPI_Isend(&feq[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],25,world,&requests[2]); + MPI_Irecv(&feq[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],15,world,&requests[3]); + if(typeLB == 2){ + MPI_Isend(&feqn[0][0][1][0],1,passzf,comm->procneigh[2][0],10,world,&requests[4]); + MPI_Irecv(&feqn[0][0][0][0],1,passzf,comm->procneigh[2][0],20,world,&requests[5]); + MPI_Isend(&feqn[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],20,world,&requests[6]); + MPI_Irecv(&feqn[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],10,world,&requests[7]); + } + update_periodic(1,subNbx-1,1,2,2,subNbz-2); + update_periodic(1,subNbx-1,subNby-2,subNby-1,2,subNbz-2); + MPI_Waitall(numrequests,requests,statuses); + + if(typeLB==1){ + update_periodic(1,subNbx-1,1,subNby-1,1,2); + update_periodic(1,subNbx-1,1,subNby-1,subNbz-2,subNbz-1); + }else if(typeLB==2){ + if(comm->myloc[2]==0){ + for(i=1; imyloc[2]==comm->procgrid[2]-1){ + for(i=1;imyloc[2]==0){ + MPI_Isend(&fnew[0][0][1][0],1,passzf,comm->procneigh[2][0],15,world,&req_send15); + MPI_Irecv(&fnew[0][0][0][0],1,passzf,comm->procneigh[2][0],25,world,&req_recv25); + } + + if(comm->myloc[2]==comm->procgrid[2]-1){ + MPI_Isend(&fnew[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],25,world,&req_send25); + MPI_Irecv(&fnew[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],15,world,&req_recv15); + } + if(comm->myloc[2]==0){ + MPI_Wait(&req_send15,&status); + MPI_Wait(&req_recv25,&status); + + for(i=1;imyloc[2]==comm->procgrid[2]-1){ + MPI_Wait(&req_send25,&status); + MPI_Wait(&req_recv15,&status); + + for(i=1;iprocneigh[0][0],15,world,&requests[0]); + MPI_Irecv(&feq[0][1][1][0],1,passxf,comm->procneigh[0][0],25,world,&requests[1]); + MPI_Isend(&feq[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],25,world,&requests[2]); + MPI_Irecv(&feq[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],15,world,&requests[3]); + if(typeLB == 2){ + MPI_Isend(&feqn[1][1][1][0],1,passxf,comm->procneigh[0][0],10,world,&requests[4]); + MPI_Irecv(&feqn[0][1][1][0],1,passxf,comm->procneigh[0][0],20,world,&requests[5]); + MPI_Isend(&feqn[subNbx-2][1][1][0],1,passxf,comm->procneigh[0][1],20,world,&requests[6]); + MPI_Irecv(&feqn[subNbx-1][1][1][0],1,passxf,comm->procneigh[0][1],10,world,&requests[7]); + } + update_periodic(2,subNbx-2,2,subNby-2,2,subNbz-2); + MPI_Waitall(numrequests,requests,statuses); + + for(i=0; iprocneigh[1][0],15,world,&requests[0]); + MPI_Irecv(&feq[0][0][1][0],1,passyf,comm->procneigh[1][0],25,world,&requests[1]); + MPI_Isend(&feq[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],25,world,&requests[2]); + MPI_Irecv(&feq[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],15,world,&requests[3]); + if(typeLB == 2){ + MPI_Isend(&feqn[0][1][1][0],1,passyf,comm->procneigh[1][0],10,world,&requests[4]); + MPI_Irecv(&feqn[0][0][1][0],1,passyf,comm->procneigh[1][0],20,world,&requests[5]); + MPI_Isend(&feqn[0][subNby-2][1][0],1,passyf,comm->procneigh[1][1],20,world,&requests[6]); + MPI_Irecv(&feqn[0][subNby-1][1][0],1,passyf,comm->procneigh[1][1],10,world,&requests[7]); + } + update_periodic(1,2,2,subNby-2,2,subNbz-2); + update_periodic(subNbx-2,subNbx-1,2,subNby-2,2,subNbz-2); + MPI_Waitall(numrequests,requests,statuses); + + for(i=0; iprocneigh[2][0],15,world,&requests[0]); + MPI_Irecv(&feq[0][0][0][0],1,passzf,comm->procneigh[2][0],25,world,&requests[1]); + MPI_Isend(&feq[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],25,world,&requests[2]); + MPI_Irecv(&feq[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],15,world,&requests[3]); + if(typeLB == 2){ + MPI_Isend(&feqn[0][0][1][0],1,passzf,comm->procneigh[2][0],10,world,&requests[4]); + MPI_Irecv(&feqn[0][0][0][0],1,passzf,comm->procneigh[2][0],20,world,&requests[5]); + MPI_Isend(&feqn[0][0][subNbz-2][0],1,passzf,comm->procneigh[2][1],20,world,&requests[6]); + MPI_Irecv(&feqn[0][0][subNbz-1][0],1,passzf,comm->procneigh[2][1],10,world,&requests[7]); + } + update_periodic(1,subNbx-1,1,2,2,subNbz-2); + update_periodic(1,subNbx-1,subNby-2,subNby-1,2,subNbz-2); + MPI_Waitall(numrequests,requests,statuses); + + update_periodic(1,subNbx-1,1,subNby-1,1,2); + update_periodic(1,subNbx-1,1,subNby-1,subNbz-2,subNbz-1); + } + +} + + diff --git a/src/USER-LB/fix_lb_fluid.h b/src/USER-LB/fix_lb_fluid.h new file mode 100755 index 0000000000..48536bdf6d --- /dev/null +++ b/src/USER-LB/fix_lb_fluid.h @@ -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 + diff --git a/src/USER-LB/fix_lb_momentum.cpp b/src/USER-LB/fix_lb_momentum.cpp new file mode 100755 index 0000000000..06da991a20 --- /dev/null +++ b/src/USER-LB/fix_lb_momentum.cpp @@ -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; ifixnfix; 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; ie; + + //Subtract vcm from the fluid. + for(int i=0; ifeqold; + 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; lparametercalc_full(); + +} diff --git a/src/USER-LB/fix_lb_momentum.h b/src/USER-LB/fix_lb_momentum.h new file mode 100755 index 0000000000..f9a55eba4a --- /dev/null +++ b/src/USER-LB/fix_lb_momentum.h @@ -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 diff --git a/src/USER-LB/fix_lb_pc.cpp b/src/USER-LB/fix_lb_pc.cpp new file mode 100755 index 0000000000..abd61a8f88 --- /dev/null +++ b/src/USER-LB/fix_lb_pc.cpp @@ -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; ifixnfix; 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; jone(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; iftm2v - 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; isublo[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; + + } + } + } diff --git a/src/USER-LB/fix_lb_pc.h b/src/USER-LB/fix_lb_pc.h new file mode 100755 index 0000000000..9287bfa74a --- /dev/null +++ b/src/USER-LB/fix_lb_pc.h @@ -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 diff --git a/src/USER-LB/fix_lb_rigid_pc_sphere.cpp b/src/USER-LB/fix_lb_rigid_pc_sphere.cpp new file mode 100755 index 0000000000..f6ae832b72 --- /dev/null +++ b/src/USER-LB/fix_lb_rigid_pc_sphere.cpp @@ -0,0 +1,1691 @@ +/* ---------------------------------------------------------------------- + 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 (version from 2008). +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "fix_lb_rigid_pc_sphere.h" +#include "atom.h" +#include "atom_vec.h" +#include "domain.h" +#include "update.h" +#include "respa.h" +#include "modify.h" +#include "group.h" +#include "comm.h" +#include "force.h" +#include "output.h" +#include "memory.h" +#include "error.h" +#include "fix_lb_fluid.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +#define MIN(A,B) ((A) < (B)) ? (A) : (B) +#define MAX(A,B) ((A) > (B)) ? (A) : (B) +/* -------------------------------------------------------------------------- */ + +FixLbRigidPCSphere::FixLbRigidPCSphere(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + int i, ibody; + + scalar_flag = 1; + extscalar = 0; + time_integrate = 1; + rigid_flag = 1; + create_attribute = 1; + virial_flag = 1; + + // perform initial allocation of atom-based arrays + // register with Atom class + + body = NULL; + up = NULL; + grow_arrays(atom->nmax); + atom->add_callback(0); + + // by default assume all of the particles interact with the fluid. + inner_nodes = 0; + + // parse command-line args + // set nbody and body[i] for each atom + + if (narg < 4) error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command"); + + // single rigid body + // nbody = 1 + // all atoms in fix group are part of body + + int iarg; + + if (strcmp(arg[3],"single") == 0) { + iarg = 4; + nbody = 1; + + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (i = 0; i < nlocal; i++) { + body[i] = -1; + if (mask[i] & groupbit) body[i] = 0; + } + + // each molecule in fix group is a rigid body + // maxmol = largest molecule # + // ncount = # of atoms in each molecule (have to sum across procs) + // nbody = # of non-zero ncount values + // use nall as incremented ptr to set body[] values for each atom + + } else if (strcmp(arg[3],"molecule") == 0) { + iarg = 4; + if (atom->molecular == 0) + error->all(FLERR,"Must use a molecular atom style with fix lb/rigid/pc/sphere molecule"); + + int *mask = atom->mask; + int *molecule = atom->molecule; + int nlocal = atom->nlocal; + + int maxmol = -1; + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) maxmol = MAX(maxmol,molecule[i]); + + int itmp; + MPI_Allreduce(&maxmol,&itmp,1,MPI_INT,MPI_MAX,world); + maxmol = itmp + 1; + + int *ncount = new int[maxmol]; + for (i = 0; i < maxmol; i++) ncount[i] = 0; + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) ncount[molecule[i]]++; + + int *nall = new int[maxmol]; + MPI_Allreduce(ncount,nall,maxmol,MPI_INT,MPI_SUM,world); + + nbody = 0; + for (i = 0; i < maxmol; i++) + if (nall[i]) nall[i] = nbody++; + else nall[i] = -1; + + for (i = 0; i < nlocal; i++) { + body[i] = -1; + if (mask[i] & groupbit) body[i] = nall[molecule[i]]; + } + + delete [] ncount; + delete [] nall; + + // each listed group is a rigid body + // check if all listed groups exist + // an atom must belong to fix group and listed group to be in rigid body + // error if atom belongs to more than 1 rigid body + + } else if (strcmp(arg[3],"group") == 0) { + if (narg < 5) error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command"); + nbody = atoi(arg[4]); + if (nbody <= 0) error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command"); + if (narg < 5+nbody) error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command"); + iarg = 5 + nbody; + + int *igroups = new int[nbody]; + for (ibody = 0; ibody < nbody; ibody++) { + igroups[ibody] = group->find(arg[5+ibody]); + if (igroups[ibody] == -1) + error->all(FLERR,"Could not find fix lb/rigid/pc/sphere group ID"); + } + + int *mask = atom->mask; + int nlocal = atom->nlocal; + + int flag = 0; + for (i = 0; i < nlocal; i++) { + body[i] = -1; + if (mask[i] & groupbit) + for (ibody = 0; ibody < nbody; ibody++) + if (mask[i] & group->bitmask[igroups[ibody]]) { + if (body[i] >= 0) flag = 1; + body[i] = ibody; + } + } + + int flagall; + MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); + if (flagall) + error->all(FLERR,"One or more atoms belong to multiple rigid bodies"); + + delete [] igroups; + + }else error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command"); + + // error check on nbody + + if (nbody == 0) error->all(FLERR,"No rigid bodies defined"); + + // create all nbody-length arrays + + memory->create(nrigid,nbody,"lb/rigid/pc/sphere:nrigid"); + memory->create(nrigid_shell,nbody,"lb/rigid/pc/sphere:nrigid_shell"); + memory->create(masstotal,nbody,"lb/rigid/pc/sphere:masstotal"); + memory->create(masstotal_shell,nbody,"lb/rigid/pc/sphere:masstotal_shell"); + memory->create(sphereradius,nbody,"lb/rigid/pc/sphere:sphereradius"); + memory->create(xcm,nbody,3,"lb/rigid/pc/sphere:xcm"); + memory->create(xcm_old,nbody,3,"lb/rigid/pc/sphere:xcm_old"); + memory->create(vcm,nbody,3,"lb/rigid/pc/sphere:vcm"); + memory->create(ucm,nbody,3,"lb/rigid/pc/sphere:ucm"); + memory->create(ucm_old,nbody,3,"lb/rigid/pc/sphere:ucm_old"); + memory->create(fcm,nbody,3,"lb/rigid/pc/sphere:fcm"); + memory->create(fcm_old,nbody,3,"lb/rigid/pc/sphere:fcm_old"); + memory->create(fcm_fluid,nbody,3,"lb/rigid/pc/sphere:fcm_fluid"); + memory->create(omega,nbody,3,"lb/rigid/pc/sphere:omega"); + memory->create(torque,nbody,3,"lb/rigid/pc/sphere:torque"); + memory->create(torque_old,nbody,3,"lb/rigid/pc/sphere:torque_old"); + memory->create(torque_fluid,nbody,3,"lb/rigid/pc/sphere:torque_fluid"); + memory->create(torque_fluid_old,nbody,3,"lb/rigid/pc/sphere:torque_fluid_old"); + memory->create(rotate,nbody,3,"lb/rigid/pc/sphere:rotate"); + memory->create(imagebody,nbody,"lb/rigid/pc/sphere:imagebody"); + memory->create(fflag,nbody,3,"lb/rigid/pc/sphere:fflag"); + memory->create(tflag,nbody,3,"lb/rigid/pc/sphere:tflag"); + + memory->create(sum,nbody,6,"lb/rigid/pc/sphere:sum"); + memory->create(all,nbody,6,"lb/rigid/pc/sphere:all"); + memory->create(remapflag,nbody,4,"lb/rigid/pc/sphere:remapflag"); + + + Gamma_MD = new double[nbody]; + + // initialize force/torque flags to default = 1.0 + + array_flag = 1; + size_array_rows = nbody; + size_array_cols = 15; + global_freq = 1; + extarray = 0; + + for (i = 0; i < nbody; i++) { + fflag[i][0] = fflag[i][1] = fflag[i][2] = 1.0; + tflag[i][0] = tflag[i][1] = tflag[i][2] = 1.0; + } + + // parse optional args that set fflag and tflag + + while (iarg < narg) { + if (strcmp(arg[iarg],"force") == 0) { + if (iarg+5 > narg) error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command"); + + int mlo,mhi; + force->bounds(arg[iarg+1],nbody,mlo,mhi); + + double xflag,yflag,zflag; + if (strcmp(arg[iarg+2],"off") == 0) xflag = 0.0; + else if (strcmp(arg[iarg+2],"on") == 0) xflag = 1.0; + else error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command"); + if (strcmp(arg[iarg+2],"off") == 0) yflag = 0.0; + else if (strcmp(arg[iarg+3],"on") == 0) yflag = 1.0; + else error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command"); + if (strcmp(arg[iarg+4],"off") == 0) zflag = 0.0; + else if (strcmp(arg[iarg+4],"on") == 0) zflag = 1.0; + else error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command"); + + int count = 0; + for (int m = mlo; m <= mhi; m++) { + fflag[m-1][0] = xflag; + fflag[m-1][1] = yflag; + fflag[m-1][2] = zflag; + count++; + } + if (count == 0) error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command"); + + iarg += 5; + } else if (strcmp(arg[iarg],"torque") == 0) { + if (iarg+5 > narg) error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command"); + + int mlo,mhi; + force->bounds(arg[iarg+1],nbody,mlo,mhi); + + double xflag,yflag,zflag; + if (strcmp(arg[iarg+2],"off") == 0) xflag = 0.0; + else if (strcmp(arg[iarg+2],"on") == 0) xflag = 1.0; + else error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command"); + if (strcmp(arg[iarg+3],"off") == 0) yflag = 0.0; + else if (strcmp(arg[iarg+3],"on") == 0) yflag = 1.0; + else error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command"); + if (strcmp(arg[iarg+4],"off") == 0) zflag = 0.0; + else if (strcmp(arg[iarg+4],"on") == 0) zflag = 1.0; + else error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command"); + + int count = 0; + for (int m = mlo; m <= mhi; m++) { + tflag[m-1][0] = xflag; + tflag[m-1][1] = yflag; + tflag[m-1][2] = zflag; + count++; + } + if (count == 0) error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command"); + + iarg += 5; + // specify if certain particles are inside the rigid spherical body, + // and therefore should not + } else if(strcmp(arg[iarg],"innerNodes")==0){ + inner_nodes = 1; + igroupinner = group->find(arg[iarg+1]); + if(igroupinner == -1) + error->all(FLERR,"Could not find fix lb/rigid/pc/sphere innerNodes group ID"); + iarg += 2; + } else error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command"); + } + + // initialize vector output quantities in case accessed before run + + for (i = 0; i < nbody; i++) { + xcm[i][0] = xcm[i][1] = xcm[i][2] = 0.0; + xcm_old[i][0] = xcm_old[i][1] = xcm_old[i][2] = 0.0; + vcm[i][0] = vcm[i][1] = vcm[i][2] = 0.0; + ucm[i][0] = ucm[i][1] = ucm[i][2] = 0.0; + ucm_old[i][0] = ucm_old[i][1] = ucm_old[i][2] = 0.0; + fcm[i][0] = fcm[i][1] = fcm[i][2] = 0.0; + fcm_old[i][0] = fcm_old[i][1] = fcm_old[i][2] = 0.0; + fcm_fluid[i][0] = fcm_fluid[i][1] = fcm_fluid[i][2] = 0.0; + torque[i][0] = torque[i][1] = torque[i][2] = 0.0; + torque_old[i][0] = torque_old[i][1] = torque_old[i][2] = 0.0; + torque_fluid[i][0] = torque_fluid[i][1] = torque_fluid[i][2] = 0.0; + torque_fluid_old[i][0] = torque_fluid_old[i][1] = torque_fluid_old[i][2] = 0.0; + } + + // nrigid[n] = # of atoms in Nth rigid body + // error if one or zero atoms + + int *ncount = new int[nbody]; + for (ibody = 0; ibody < nbody; ibody++) ncount[ibody] = 0; + + int nlocal = atom->nlocal; + + for (i = 0; i < nlocal; i++) + if (body[i] >= 0) ncount[body[i]]++; + + MPI_Allreduce(ncount,nrigid,nbody,MPI_INT,MPI_SUM,world); + + //count the number of atoms in the shell. + if(inner_nodes == 1){ + int *mask = atom->mask; + for(ibody=0; ibodybitmask[igroupinner])){ + if(body[i] >= 0) ncount[body[i]]++; + } + } + + MPI_Allreduce(ncount,nrigid_shell,nbody,MPI_INT,MPI_SUM,world); + }else { + for(ibody=0; ibody < nbody; ibody++) nrigid_shell[ibody]=nrigid[ibody]; + } + + delete [] ncount; + + for (ibody = 0; ibody < nbody; ibody++) + if (nrigid[ibody] <= 1) error->all(FLERR,"One or zero atoms in rigid body"); + + // set image flags for each rigid body to default values + // will be reset during init() based on xcm and then by pre_neighbor() + // set here, so image value will persist from run to run + + for (ibody = 0; ibody < nbody; ibody++) + imagebody[ibody] = ((tagint) IMGMAX << IMG2BITS) | + ((tagint) IMGMAX << IMGBITS) | IMGMAX; + + // print statistics + + int nsum = 0; + for (ibody = 0; ibody < nbody; ibody++) nsum += nrigid[ibody]; + + if (comm->me == 0) { + if (screen) fprintf(screen,"%d rigid bodies with %d atoms\n",nbody,nsum); + if (logfile) fprintf(logfile,"%d rigid bodies with %d atoms\n",nbody,nsum); + } + + int groupbit_lb_fluid = 0; + + for(int ifix=0; ifixnfix; 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/rigid/pc/sphere fix"); + + int *mask = atom->mask; + if(inner_nodes == 1){ + for(int j=0; jbitmask[igroupinner]) && !(mask[j] & groupbit_lb_fluid)) + error->one(FLERR,"use the innerNodes keyword in the lb/rigid/pc/sphere fix for atoms which do not interact with the lb/fluid"); + + // If inner nodes are present, which should not interact with the fluid, make + // sure these are not used by the lb/fluid fix to apply a force to the fluid. + if((mask[j] & groupbit) && (mask[j] & groupbit_lb_fluid) && (mask[j] & group->bitmask[igroupinner])) + error->one(FLERR,"the inner nodes specified in lb/rigid/pc/sphere should not be included in the lb/fluid fix"); + } + }else{ + for(int j=0; jone(FLERR,"use the innerNodes keyword in the lb/rigid/pc/sphere fix for atoms which do not interact with the lb/fluid"); + } + } + +} + +/* ---------------------------------------------------------------------- */ + +FixLbRigidPCSphere::~FixLbRigidPCSphere() +{ + // unregister callbacks to this fix from Atom class + + atom->delete_callback(id,0); + + // delete locally stored arrays + + memory->destroy(body); + memory->destroy(up); + + // delete nbody-length arrays + + memory->destroy(nrigid); + memory->destroy(nrigid_shell); + memory->destroy(masstotal); + memory->destroy(masstotal_shell); + memory->destroy(sphereradius); + memory->destroy(xcm); + memory->destroy(xcm_old); + memory->destroy(vcm); + memory->destroy(ucm); + memory->destroy(ucm_old); + memory->destroy(fcm); + memory->destroy(fcm_old); + memory->destroy(fcm_fluid); + memory->destroy(omega); + memory->destroy(torque); + memory->destroy(torque_old); + memory->destroy(torque_fluid); + memory->destroy(torque_fluid_old); + memory->destroy(rotate); + memory->destroy(imagebody); + memory->destroy(fflag); + memory->destroy(tflag); + + memory->destroy(sum); + memory->destroy(all); + memory->destroy(remapflag); + + delete [] Gamma_MD; +} + +/* ---------------------------------------------------------------------- */ + +int FixLbRigidPCSphere::setmask() +{ + int mask = 0; + mask |= INITIAL_INTEGRATE; + mask |= FINAL_INTEGRATE; + mask |= PRE_NEIGHBOR; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixLbRigidPCSphere::init() +{ + int i,itype,ibody; + + // warn if more than one rigid fix + + int count = 0; + for (int i = 0; i < modify->nfix; i++) + if (strcmp(modify->fix[i]->style,"lb/rigid/pc/sphere") == 0) count++; + if (count > 1 && comm->me == 0) error->warning(FLERR,"More than one fix lb/rigid/pc/sphere"); + + // timestep info + + dtv = update->dt; + dtf = 0.5 * update->dt * force->ftm2v; + + int *type = atom->type; + int nlocal = atom->nlocal; + double **x = atom->x; + tagint *image = atom->image; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *periodicity = domain->periodicity; + double xprd = domain->xprd; + double yprd = domain->yprd; + double zprd = domain->zprd; + + double **mu = atom->mu; + double *radius = atom->radius; + int *ellipsoid = atom->ellipsoid; + int extended = 0; + int *mask = atom->mask; + + // Warn if any extended particles are included. + if (atom->radius_flag || atom->ellipsoid_flag || atom->mu_flag) { + int flag = 0; + for (i = 0; i < nlocal; i++) { + if (body[i] < 0) continue; + if (radius && radius[i] > 0.0) flag = 1; + if (ellipsoid && ellipsoid[i] >= 0) flag = 1; + if (mu && mu[i][3] > 0.0) flag = 1; + } + + MPI_Allreduce(&flag,&extended,1,MPI_INT,MPI_MAX,world); + } + if(extended) + error->warning(FLERR,"Fix lb/rigid/pc/sphere assumes point particles"); + + // compute masstotal & center-of-mass of each rigid body + + for (ibody = 0; ibody < nbody; ibody++) + for (i = 0; i < 6; i++) sum[ibody][i] = 0.0; + int xbox,ybox,zbox; + double massone,xunwrap,yunwrap,zunwrap; + + for (i = 0; i < nlocal; i++) { + if (body[i] < 0) continue; + ibody = body[i]; + + xbox = (image[i] & IMGMASK) - IMGMAX; + ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; + zbox = (image[i] >> IMG2BITS) - IMGMAX; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + + if ((xbox && !periodicity[0]) || (ybox && !periodicity[1]) || + (zbox && !periodicity[2])) + error->one(FLERR,"Fix lb/rigid/pc/sphere atom has non-zero image flag " + "in a non-periodic dimension"); + + xunwrap = x[i][0] + xbox*xprd; + yunwrap = x[i][1] + ybox*yprd; + zunwrap = x[i][2] + zbox*zprd; + + sum[ibody][0] += xunwrap * massone; + sum[ibody][1] += yunwrap * massone; + sum[ibody][2] += zunwrap * massone; + sum[ibody][3] += massone; + if(inner_nodes == 1){ + if(!(mask[i] & group->bitmask[igroupinner])){ + sum[ibody][4] += massone; + } + }else{ + sum[ibody][4] += massone; + } + } + + MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world); + + for (ibody = 0; ibody < nbody; ibody++) { + masstotal[ibody] = all[ibody][3]; + masstotal_shell[ibody] = all[ibody][4]; + xcm[ibody][0] = all[ibody][0]/masstotal[ibody]; + xcm[ibody][1] = all[ibody][1]/masstotal[ibody]; + xcm[ibody][2] = all[ibody][2]/masstotal[ibody]; + } + + // Calculate the radius of the rigid body, and assign the value for gamma: + double dx,dy,dz; + double *Gamma = fix_lb_fluid->Gamma; + double dm_lb = fix_lb_fluid->dm_lb; + double dt_lb = fix_lb_fluid->dt_lb; + + for (ibody = 0; ibody < nbody; ibody++) + for (i = 0; i < 6; i++) sum[ibody][i] = 0.0; + + for (i=0; ibitmask[igroupinner])){ + ibody = body[i]; + + xbox = (image[i] & IMGMASK) - IMGMAX; + ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; + zbox = (image[i] >> IMG2BITS) - IMGMAX; + + xunwrap = x[i][0] + xbox*xprd; + yunwrap = x[i][1] + ybox*yprd; + zunwrap = x[i][2] + zbox*zprd; + + dx = xunwrap - xcm[ibody][0]; + dy = yunwrap - xcm[ibody][1]; + dz = zunwrap - xcm[ibody][2]; + + sum[ibody][0] += dx*dx + dy*dy + dz*dz; + sum[ibody][1] += Gamma[type[i]]; + } + }else{ + ibody = body[i]; + + xbox = (image[i] & IMGMASK) - IMGMAX; + ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; + zbox = (image[i] >> IMG2BITS) - IMGMAX; + + xunwrap = x[i][0] + xbox*xprd; + yunwrap = x[i][1] + ybox*yprd; + zunwrap = x[i][2] + zbox*zprd; + + dx = xunwrap - xcm[ibody][0]; + dy = yunwrap - xcm[ibody][1]; + dz = zunwrap - xcm[ibody][2]; + + sum[ibody][0] += dx*dx + dy*dy + dz*dz; + sum[ibody][1] += Gamma[type[i]]; + } + } + + MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world); + + for(ibody=0; ibody < nbody; ibody++){ + sphereradius[ibody] = sqrt(all[ibody][0]/nrigid_shell[ibody]); + Gamma_MD[ibody] = all[ibody][1]*dm_lb/dt_lb/nrigid_shell[ibody]; + } + + // Check that all atoms in the rigid body have the same value of gamma. + double eps = 1.0e-7; + for (i=0; ibitmask[igroupinner])){ + ibody = body[i]; + + if(Gamma_MD[ibody]*dt_lb/dm_lb - Gamma[type[i]] > eps) + error->one(FLERR,"All atoms in a rigid body must have the same gamma value"); + } + }else{ + ibody = body[i]; + + if(Gamma_MD[ibody]*dt_lb/dm_lb - Gamma[type[i]] > eps) + error->one(FLERR,"All atoms in a rigid body must have the same gamma value"); + } + } + + + // remap the xcm of each body back into simulation box if needed + // only really necessary the 1st time a run is performed + + pre_neighbor(); + + + // temperature scale factor + + double ndof = 0.0; + for (ibody = 0; ibody < nbody; ibody++) { + ndof += fflag[ibody][0] + fflag[ibody][1] + fflag[ibody][2]; + ndof += tflag[ibody][0] + tflag[ibody][1] + tflag[ibody][2]; + } + if (ndof > 0.0) tfactor = force->mvv2e / (ndof * force->boltz); + else tfactor = 0.0; + +} +/* ---------------------------------------------------------------------- */ + +void FixLbRigidPCSphere::setup(int vflag) +{ + int i,n,ibody; + double massone,radone; + + // vcm = velocity of center-of-mass of each rigid body + // fcm = force on center-of-mass of each rigid body + + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; + int nlocal = atom->nlocal; + + tagint *image = atom->image; + int *periodicity = domain->periodicity; + + double unwrap[3]; + double dx,dy,dz; + + + + for (ibody = 0; ibody < nbody; ibody++) + for (i = 0; i < 6; i++) sum[ibody][i] = 0.0; + + for (i = 0; i < nlocal; i++) { + if (body[i] < 0) continue; + ibody = body[i]; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + + sum[ibody][0] += v[i][0] * massone; + sum[ibody][1] += v[i][1] * massone; + sum[ibody][2] += v[i][2] * massone; + sum[ibody][3] += f[i][0]; + sum[ibody][4] += f[i][1]; + sum[ibody][5] += f[i][2]; + } + + MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world); + + for (ibody = 0; ibody < nbody; ibody++) { + vcm[ibody][0] = all[ibody][0]/masstotal[ibody]; + vcm[ibody][1] = all[ibody][1]/masstotal[ibody]; + vcm[ibody][2] = all[ibody][2]/masstotal[ibody]; + fcm[ibody][0] = all[ibody][3]; + fcm[ibody][1] = all[ibody][4]; + fcm[ibody][2] = all[ibody][5]; + } + + // omega = angular velocity of each rigid body + // Calculated as the average of the angular velocities of the + // individual atoms comprising the rigid body. + // torque = torque on each rigid body + for (ibody = 0; ibody < nbody; ibody++) + for (i = 0; i < 6; i++) sum[ibody][i] = 0.0; + + for (i = 0; i < nlocal; i++) { + if (body[i] < 0) continue; + ibody = body[i]; + + domain->unmap(x[i],image[i],unwrap); + + dx = unwrap[0] - xcm[ibody][0]; + dy = unwrap[1] - xcm[ibody][1]; + dz = unwrap[2] - xcm[ibody][2]; + + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + + sum[ibody][0] += (dy * (v[i][2]-vcm[ibody][2]) - dz * (v[i][1]-vcm[ibody][1]))/(dx*dx+dy*dy+dz*dz); + sum[ibody][1] += (dz * (v[i][0]-vcm[ibody][0]) - dx * (v[i][2]-vcm[ibody][2]))/(dx*dx+dy*dy+dz*dz); + sum[ibody][2] += (dx * (v[i][1]-vcm[ibody][1]) - dy * (v[i][0]-vcm[ibody][0]))/(dx*dx+dy*dy+dz*dz); + sum[ibody][3] += dy * f[i][2] - dz * f[i][1]; + sum[ibody][4] += dz * f[i][0] - dx * f[i][2]; + sum[ibody][5] += dx * f[i][1] - dy * f[i][0]; + } + + MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world); + + for (ibody = 0; ibody < nbody; ibody++) { + omega[ibody][0] = all[ibody][0]/nrigid[ibody]; + omega[ibody][1] = all[ibody][1]/nrigid[ibody]; + omega[ibody][2] = all[ibody][2]/nrigid[ibody]; + torque[ibody][0] = all[ibody][3]; + torque[ibody][1] = all[ibody][4]; + torque[ibody][2] = all[ibody][5]; + } + + // virial setup before call to set_v + + if (vflag) v_setup(vflag); + else evflag = 0; + + // Set the velocities + set_v(); + + if (vflag_global) + for (n = 0; n < 6; n++) virial[n] *= 2.0; + if (vflag_atom) { + for (i = 0; i < nlocal; i++) + for (n = 0; n < 6; n++) + vatom[i][n] *= 2.0; + } + +} + +/* ---------------------------------------------------------------------- */ + +void FixLbRigidPCSphere::initial_integrate(int vflag) +{ + double dtfm; + + int i,n,ibody; + + double massone,radone; + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; + int nlocal = atom->nlocal; + + tagint *image = atom->image; + + double unwrap[3]; + double dx,dy,dz; + + int *mask = atom->mask; + + // compute the fluid velocity at the initial particle positions + compute_up(); + + for (ibody = 0; ibody < nbody; ibody++) + for (i = 0; i < 6; i++) sum[ibody][i] = 0.0; + // Store the fluid velocity at the center of mass + + for (i = 0; i < nlocal; i++) { + if (body[i] < 0) continue; + ibody = body[i]; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + + if(inner_nodes == 1){ + if(!(mask[i] & group->bitmask[igroupinner])){ + sum[ibody][0] += up[i][0]*massone; + sum[ibody][1] += up[i][1]*massone; + sum[ibody][2] += up[i][2]*massone; + } + }else{ + sum[ibody][0] += up[i][0]*massone; + sum[ibody][1] += up[i][1]*massone; + sum[ibody][2] += up[i][2]*massone; + } + } + MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world); + + for (ibody = 0; ibody < nbody; ibody++) { + ucm[ibody][0] = all[ibody][0]/masstotal_shell[ibody]; + ucm[ibody][1] = all[ibody][1]/masstotal_shell[ibody]; + ucm[ibody][2] = all[ibody][2]/masstotal_shell[ibody]; + } + + //Store the total torque due to the fluid. + for (ibody = 0; ibody < nbody; ibody++) + for(i = 0; i < 6; i++) sum[ibody][i] = 0.0; + + for(i = 0; iunmap(x[i],image[i],unwrap); + + dx = unwrap[0] - xcm[ibody][0]; + dy = unwrap[1] - xcm[ibody][1]; + dz = unwrap[2] - xcm[ibody][2]; + + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + + if(inner_nodes == 1){ + if(!(mask[i] & group->bitmask[igroupinner])){ + sum[ibody][0] += Gamma_MD[ibody]*(dy * ((up[i][2]-vcm[ibody][2])) - + dz * ((up[i][1]-vcm[ibody][1]))); + sum[ibody][1] += Gamma_MD[ibody]*(dz * ((up[i][0]-vcm[ibody][0])) - + dx * ((up[i][2]-vcm[ibody][2]))); + sum[ibody][2] += Gamma_MD[ibody]*(dx * ((up[i][1]-vcm[ibody][1])) - + dy * ((up[i][0]-vcm[ibody][0]))); + sum[ibody][3] += -Gamma_MD[ibody]*(v[i][0]-up[i][0]); + sum[ibody][4] += -Gamma_MD[ibody]*(v[i][1]-up[i][1]); + sum[ibody][5] += -Gamma_MD[ibody]*(v[i][2]-up[i][2]); + } + }else{ + sum[ibody][0] += Gamma_MD[ibody]*(dy * ((up[i][2]-vcm[ibody][2])) - + dz * ((up[i][1]-vcm[ibody][1]))); + sum[ibody][1] += Gamma_MD[ibody]*(dz * ((up[i][0]-vcm[ibody][0])) - + dx * ((up[i][2]-vcm[ibody][2]))); + sum[ibody][2] += Gamma_MD[ibody]*(dx * ((up[i][1]-vcm[ibody][1])) - + dy * ((up[i][0]-vcm[ibody][0]))); + sum[ibody][3] += -Gamma_MD[ibody]*(v[i][0]-up[i][0]); + sum[ibody][4] += -Gamma_MD[ibody]*(v[i][1]-up[i][1]); + sum[ibody][5] += -Gamma_MD[ibody]*(v[i][2]-up[i][2]); + } + } + + MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world); + + for (ibody = 0; ibody < nbody; ibody++) { + torque_fluid[ibody][0] = all[ibody][0]; + torque_fluid[ibody][1] = all[ibody][1]; + torque_fluid[ibody][2] = all[ibody][2]; + fcm_fluid[ibody][0] = all[ibody][3]; + fcm_fluid[ibody][1] = all[ibody][4]; + fcm_fluid[ibody][2] = all[ibody][5]; + } + + for (int ibody = 0; ibody < nbody; ibody++) { + fcm_old[ibody][0] = fcm[ibody][0]; + fcm_old[ibody][1] = fcm[ibody][1]; + fcm_old[ibody][2] = fcm[ibody][2]; + torque_old[ibody][0] = torque[ibody][0]; + torque_old[ibody][1] = torque[ibody][1]; + torque_old[ibody][2] = torque[ibody][2]; + torque_fluid_old[ibody][0] = torque_fluid[ibody][0]; + torque_fluid_old[ibody][1] = torque_fluid[ibody][1]; + torque_fluid_old[ibody][2] = torque_fluid[ibody][2]; + ucm_old[ibody][0] = ucm[ibody][0]; + ucm_old[ibody][1] = ucm[ibody][1]; + ucm_old[ibody][2] = ucm[ibody][2]; + xcm_old[ibody][0] = xcm[ibody][0]; + xcm_old[ibody][1] = xcm[ibody][1]; + xcm_old[ibody][2] = xcm[ibody][2]; + + // update xcm by full step + + dtfm = dtf / masstotal[ibody]; + xcm[ibody][0] += dtv * vcm[ibody][0]+(fcm[ibody][0]+fcm_fluid[ibody][0]/force->ftm2v)*fflag[ibody][0]*dtfm*dtv; + xcm[ibody][1] += dtv * vcm[ibody][1]+(fcm[ibody][1]+fcm_fluid[ibody][1]/force->ftm2v)*fflag[ibody][1]*dtfm*dtv; + xcm[ibody][2] += dtv * vcm[ibody][2]+(fcm[ibody][2]+fcm_fluid[ibody][2]/force->ftm2v)*fflag[ibody][2]*dtfm*dtv; + + rotate[ibody][0] = omega[ibody][0]*dtv + tflag[ibody][0]*(torque[ibody][0]*force->ftm2v+torque_fluid[ibody][0])* + dtv*dtv*5.0/(4.0*masstotal[ibody]*sphereradius[ibody]*sphereradius[ibody]); + rotate[ibody][1] = omega[ibody][1]*dtv + tflag[ibody][1]*(torque[ibody][1]*force->ftm2v+torque_fluid[ibody][1])* + dtv*dtv*5.0/(4.0*masstotal[ibody]*sphereradius[ibody]*sphereradius[ibody]); + rotate[ibody][2] = omega[ibody][2]*dtv + tflag[ibody][2]*(torque[ibody][2]*force->ftm2v+torque_fluid[ibody][2])* + dtv*dtv*5.0/(4.0*masstotal[ibody]*sphereradius[ibody]*sphereradius[ibody]); + + // Approximate vcm + expminusdttimesgamma = exp(-Gamma_MD[ibody]*dtv*nrigid_shell[ibody]/masstotal[ibody]); + force_factor = force->ftm2v/Gamma_MD[ibody]/nrigid_shell[ibody]; + + if(fflag[ibody][0]==1){ + vcm[ibody][0] = expminusdttimesgamma*(vcm[ibody][0] - ucm[ibody][0] - fcm[ibody][0]*force_factor) + + ucm[ibody][0] + fcm[ibody][0]*force_factor; + } + if(fflag[ibody][1]==1){ + vcm[ibody][1] = expminusdttimesgamma*(vcm[ibody][1] - ucm[ibody][1] - fcm[ibody][1]*force_factor) + + ucm[ibody][1] + fcm[ibody][1]*force_factor; + } + if(fflag[ibody][2]==1){ + vcm[ibody][2] = expminusdttimesgamma*(vcm[ibody][2] - ucm[ibody][2] - fcm[ibody][2]*force_factor) + + ucm[ibody][2] + fcm[ibody][2]*force_factor; + } + + // Approximate angmom + torque_factor = 5.0*Gamma_MD[ibody]*nrigid_shell[ibody]/(3.0*masstotal[ibody]); + expminusdttimesgamma = exp(-dtv*torque_factor); + + + if(tflag[ibody][0]==1){ + omega[ibody][0] = expminusdttimesgamma*(omega[ibody][0] - (3.0/(2.0*nrigid_shell[ibody]*sphereradius[ibody]*sphereradius[ibody]*Gamma_MD[ibody]))* + (force->ftm2v*torque[ibody][0] + torque_fluid[ibody][0])) + + (3.0/(2.0*nrigid_shell[ibody]*sphereradius[ibody]*sphereradius[ibody]*Gamma_MD[ibody]))* + (force->ftm2v*torque[ibody][0] + torque_fluid[ibody][0]); + } + if(tflag[ibody][1]==1){ + omega[ibody][1] = expminusdttimesgamma*(omega[ibody][1] - (3.0/(2.0*nrigid_shell[ibody]*sphereradius[ibody]*sphereradius[ibody]*Gamma_MD[ibody]))* + (force->ftm2v*torque[ibody][1] + torque_fluid[ibody][1])) + + (3.0/(2.0*nrigid_shell[ibody]*sphereradius[ibody]*sphereradius[ibody]*Gamma_MD[ibody]))* + (force->ftm2v*torque[ibody][1] + torque_fluid[ibody][1]); + } + if(tflag[ibody][2]==1){ + omega[ibody][2] = expminusdttimesgamma*(omega[ibody][2] - (3.0/(2.0*nrigid_shell[ibody]*sphereradius[ibody]*sphereradius[ibody]*Gamma_MD[ibody]))* + (force->ftm2v*torque[ibody][2] + torque_fluid[ibody][2])) + + (3.0/(2.0*nrigid_shell[ibody]*sphereradius[ibody]*sphereradius[ibody]*Gamma_MD[ibody]))* + (force->ftm2v*torque[ibody][2] + torque_fluid[ibody][2]); + } + + } + // virial setup before call to set_xv + + if (vflag) v_setup(vflag); + else evflag = 0; + + set_xv(); + +} + + +/* ---------------------------------------------------------------------- */ + +void FixLbRigidPCSphere::final_integrate() +{ + int i,ibody; + double xy,xz,yz; + + // sum over atoms to get force and torque on rigid body + double massone; + tagint *image = atom->image; + double **x = atom->x; + double **f = atom->f; + double **v = atom->v; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; + int nlocal = atom->nlocal; + + double unwrap[3]; + double dx,dy,dz; + + int *mask = atom->mask; + + for (ibody = 0; ibody < nbody; ibody++) + for (i = 0; i < 6; i++) sum[ibody][i] = 0.0; + + for (i = 0; i < nlocal; i++) { + if (body[i] < 0) continue; + ibody = body[i]; + + sum[ibody][0] += f[i][0]; + sum[ibody][1] += f[i][1]; + sum[ibody][2] += f[i][2]; + + domain->unmap(x[i],image[i],unwrap); + + dx = unwrap[0] - xcm[ibody][0]; + dy = unwrap[1] - xcm[ibody][1]; + dz = unwrap[2] - xcm[ibody][2]; + + sum[ibody][3] += dy*f[i][2] - dz*f[i][1]; + sum[ibody][4] += dz*f[i][0] - dx*f[i][2]; + sum[ibody][5] += dx*f[i][1] - dy*f[i][0]; + } + + MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world); + + //Compute the correction to the velocity and angular momentum due to the non-fluid forces: + for (ibody = 0; ibody < nbody; ibody++) { + fcm[ibody][0] = all[ibody][0]; + fcm[ibody][1] = all[ibody][1]; + fcm[ibody][2] = all[ibody][2]; + torque[ibody][0] = all[ibody][3]; + torque[ibody][1] = all[ibody][4]; + torque[ibody][2] = all[ibody][5]; + + expminusdttimesgamma = exp(-dtv*Gamma_MD[ibody]*nrigid_shell[ibody]/masstotal[ibody]); + DMDcoeff= (dtv - (masstotal[ibody]/nrigid_shell[ibody])*(1.0-expminusdttimesgamma)/Gamma_MD[ibody]); + + vcm[ibody][0] += fflag[ibody][0]*DMDcoeff*force->ftm2v*(fcm[ibody][0]-fcm_old[ibody][0])/Gamma_MD[ibody]/dtv/nrigid_shell[ibody]; + vcm[ibody][1] += fflag[ibody][1]*DMDcoeff*force->ftm2v*(fcm[ibody][1]-fcm_old[ibody][1])/Gamma_MD[ibody]/dtv/nrigid_shell[ibody]; + vcm[ibody][2] += fflag[ibody][2]*DMDcoeff*force->ftm2v*(fcm[ibody][2]-fcm_old[ibody][2])/Gamma_MD[ibody]/dtv/nrigid_shell[ibody]; + + torque_factor = 5.0*Gamma_MD[ibody]*nrigid_shell[ibody]/(3.0*masstotal[ibody]); + expminusdttimesgamma = exp(-dtv*torque_factor); + DMDcoeff = (dtv - (1.0-expminusdttimesgamma)/torque_factor); + + omega[ibody][0] += tflag[ibody][0]*(3.0/(2.0*nrigid_shell[ibody]*sphereradius[ibody]*sphereradius[ibody]*Gamma_MD[ibody]))*DMDcoeff* + force->ftm2v*(torque[ibody][0] - torque_old[ibody][0])/dtv; + omega[ibody][1] += tflag[ibody][1]*(3.0/(2.0*nrigid_shell[ibody]*sphereradius[ibody]*sphereradius[ibody]*Gamma_MD[ibody]))*DMDcoeff* + force->ftm2v*(torque[ibody][1] - torque_old[ibody][1])/dtv; + omega[ibody][2] += tflag[ibody][2]*(3.0/(2.0*nrigid_shell[ibody]*sphereradius[ibody]*sphereradius[ibody]*Gamma_MD[ibody]))*DMDcoeff* + force->ftm2v*(torque[ibody][2] - torque_old[ibody][2])/dtv; + } + + //Next, calculate the correction to the velocity and angular momentum due to the fluid forces: + //Calculate the fluid velocity at the new particle locations. + compute_up(); + + // store fluid quantities for the total body + for (ibody = 0; ibody < nbody; ibody++) + for (i = 0; i < 6; i++) sum[ibody][i] = 0.0; + // Store the fluid velocity at the center of mass, and the total force + // due to the fluid. + for (i = 0; i < nlocal; i++) { + if (body[i] < 0) continue; + ibody = body[i]; + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + + domain->unmap(x[i],image[i],unwrap); + + dx = unwrap[0] - xcm[ibody][0]; + dy = unwrap[1] - xcm[ibody][1]; + dz = unwrap[2] - xcm[ibody][2]; + + if(inner_nodes == 1){ + if(!(mask[i] & group->bitmask[igroupinner])){ + sum[ibody][0] += up[i][0]*massone; + sum[ibody][1] += up[i][1]*massone; + sum[ibody][2] += up[i][2]*massone; + sum[ibody][3] += Gamma_MD[ibody]*(dy * ((up[i][2]-vcm[ibody][2])) - + dz * ((up[i][1]-vcm[ibody][1]))); + sum[ibody][4] += Gamma_MD[ibody]*(dz * ((up[i][0]-vcm[ibody][0])) - + dx * ((up[i][2]-vcm[ibody][2]))); + sum[ibody][5] += Gamma_MD[ibody]*(dx * ((up[i][1]-vcm[ibody][1])) - + dy * ((up[i][0]-vcm[ibody][0]))); + } + }else{ + sum[ibody][0] += up[i][0]*massone; + sum[ibody][1] += up[i][1]*massone; + sum[ibody][2] += up[i][2]*massone; + sum[ibody][3] += Gamma_MD[ibody]*(dy * ((up[i][2]-vcm[ibody][2])) - + dz * ((up[i][1]-vcm[ibody][1]))); + sum[ibody][4] += Gamma_MD[ibody]*(dz * ((up[i][0]-vcm[ibody][0])) - + dx * ((up[i][2]-vcm[ibody][2]))); + sum[ibody][5] += Gamma_MD[ibody]*(dx * ((up[i][1]-vcm[ibody][1])) - + dy * ((up[i][0]-vcm[ibody][0]))); + } + } + + MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world); + + for (ibody = 0; ibody < nbody; ibody++) { + ucm[ibody][0] = all[ibody][0]/masstotal_shell[ibody]; + ucm[ibody][1] = all[ibody][1]/masstotal_shell[ibody]; + ucm[ibody][2] = all[ibody][2]/masstotal_shell[ibody]; + torque_fluid[ibody][0] = all[ibody][3]; + torque_fluid[ibody][1] = all[ibody][4]; + torque_fluid[ibody][2] = all[ibody][5]; + } + + for (ibody = 0; ibody < nbody; ibody++) { + + expminusdttimesgamma = exp(-dtv*Gamma_MD[ibody]*nrigid_shell[ibody]/masstotal[ibody]); + DMDcoeff= (dtv - (masstotal[ibody]/nrigid_shell[ibody])*(1.0-expminusdttimesgamma)/Gamma_MD[ibody]); + + vcm[ibody][0] += DMDcoeff*fflag[ibody][0]*(ucm[ibody][0]-ucm_old[ibody][0])/dtv; + vcm[ibody][1] += DMDcoeff*fflag[ibody][1]*(ucm[ibody][1]-ucm_old[ibody][1])/dtv; + vcm[ibody][2] += DMDcoeff*fflag[ibody][2]*(ucm[ibody][2]-ucm_old[ibody][2])/dtv; + + torque_factor = 5.0*Gamma_MD[ibody]*nrigid_shell[ibody]/(3.0*masstotal[ibody]); + expminusdttimesgamma = exp(-dtv*torque_factor); + DMDcoeff = (dtv - (1.0-expminusdttimesgamma)/torque_factor); + + omega[ibody][0] += tflag[ibody][0]*(3.0/(2.0*nrigid_shell[ibody]*sphereradius[ibody]*sphereradius[ibody]*Gamma_MD[ibody]))* + DMDcoeff*(torque_fluid[ibody][0] - torque_fluid_old[ibody][0])/dtv; + omega[ibody][1] += tflag[ibody][1]*(3.0/(2.0*nrigid_shell[ibody]*sphereradius[ibody]*sphereradius[ibody]*Gamma_MD[ibody]))* + DMDcoeff*(torque_fluid[ibody][1] - torque_fluid_old[ibody][1])/dtv; + omega[ibody][2] += tflag[ibody][2]*(3.0/(2.0*nrigid_shell[ibody]*sphereradius[ibody]*sphereradius[ibody]*Gamma_MD[ibody]))* + DMDcoeff*(torque_fluid[ibody][2] - torque_fluid_old[ibody][2])/dtv; + } + + set_v(); + +} + +/* ---------------------------------------------------------------------- + set space-frame velocity of each atom in a rigid body + set omega and angmom of extended particles + v = Vcm + (W cross (x - Xcm)) +------------------------------------------------------------------------- */ + +void FixLbRigidPCSphere::set_v() +{ + int ibody; + int xbox,ybox,zbox; + double dx,dy,dz; + double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone; + double vr[6]; + + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; + tagint *image = atom->image; + int nlocal = atom->nlocal; + + double xprd = domain->xprd; + double yprd = domain->yprd; + double zprd = domain->zprd; + double xunwrap,yunwrap,zunwrap; + + // set v of each atom + + for (int i = 0; i < nlocal; i++) { + if (body[i] < 0) continue; + ibody = body[i]; + + xbox = (image[i] & IMGMASK) - IMGMAX; + ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; + zbox = (image[i] >> IMG2BITS) - IMGMAX; + + xunwrap = x[i][0] + xbox*xprd; + yunwrap = x[i][1] + ybox*yprd; + zunwrap = x[i][2] + zbox*zprd; + + dx = xunwrap - xcm[ibody][0]; + dy = yunwrap - xcm[ibody][1]; + dz = zunwrap - xcm[ibody][2]; + + // save old velocities for virial. + if(evflag){ + v0 = v[i][0]; + v1 = v[i][1]; + v2 = v[i][2]; + } + + v[i][0] = (omega[ibody][1]*dz - omega[ibody][2]*dy) + vcm[ibody][0]; + v[i][1] = (omega[ibody][2]*dx - omega[ibody][0]*dz) + vcm[ibody][1]; + v[i][2] = (omega[ibody][0]*dy - omega[ibody][1]*dx) + vcm[ibody][2]; + + // virial = unwrapped coords dotted into body constraint force + // body constraint force = implied force due to v change minus f external + // assume f does not include forces internal to body + // 1/2 factor b/c initial_integrate contributes other half + // assume per-atom contribution is due to constraint force on that atom + + if (evflag) { + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + fc0 = massone*(v[i][0] - v0)/dtf - f[i][0] + Gamma_MD[ibody]*(v0-up[i][0]); + fc1 = massone*(v[i][1] - v1)/dtf - f[i][1] + Gamma_MD[ibody]*(v1-up[i][1]); + fc2 = massone*(v[i][2] - v2)/dtf - f[i][2] + Gamma_MD[ibody]*(v2-up[i][2]); + + xbox = (image[i] & IMGMASK) - IMGMAX; + ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; + zbox = (image[i] >> IMG2BITS) - IMGMAX; + + x0 = x[i][0] + xbox*xprd; + x1 = x[i][1] + ybox*yprd; + x2 = x[i][2] + zbox*zprd; + + vr[0] = 0.5*x0*fc0; + vr[1] = 0.5*x1*fc1; + vr[2] = 0.5*x2*fc2; + vr[3] = 0.5*x0*fc1; + vr[4] = 0.5*x0*fc2; + vr[5] = 0.5*x1*fc2; + + v_tally(1,&i,1.0,vr); + } + + + } + +} +/* ---------------------------------------------------------------------- + set space-frame coords and velocity of each atom in each rigid body + set orientation and rotation of extended particles + x = Q displace + Xcm, mapped back to periodic box + v = Vcm + (W cross (x - Xcm)) +------------------------------------------------------------------------- */ + +void FixLbRigidPCSphere::set_xv() +{ + int ibody; + int xbox,ybox,zbox; + double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone; + double vr[6]; + + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; + tagint *image = atom->image; + int nlocal = atom->nlocal; + + double xprd = domain->xprd; + double yprd = domain->yprd; + double zprd = domain->zprd; + double xunwrap,yunwrap,zunwrap; + double dx,dy,dz; + + + // set x and v of each atom + + for (int i = 0; i < nlocal; i++) { + if (body[i] < 0) continue; + ibody = body[i]; + + xbox = (image[i] & IMGMASK) - IMGMAX; + ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; + zbox = (image[i] >> IMG2BITS) - IMGMAX; + + xunwrap = x[i][0] + xbox*xprd; + yunwrap = x[i][1] + ybox*yprd; + zunwrap = x[i][2] + zbox*zprd; + + dx = xunwrap - xcm_old[ibody][0]; + dy = yunwrap - xcm_old[ibody][1]; + dz = zunwrap - xcm_old[ibody][2]; + + // save old positions and velocities for virial + if(evflag){ + x0 = xunwrap; + x1 = yunwrap; + x2 = zunwrap; + v0 = v[i][0]; + v1 = v[i][1]; + v2 = v[i][2]; + } + + // x = displacement from center-of-mass, based on body orientation + // v = vcm + omega around center-of-mass + + x[i][0] = dx; + x[i][1] = dy; + x[i][2] = dz; + + //Perform the rotations: + dx = x[i][0]; + dy = x[i][1]; + dz = x[i][2]; + x[i][0] = cos(rotate[ibody][2])*dx - sin(rotate[ibody][2])*dy; + x[i][1] = sin(rotate[ibody][2])*dx + cos(rotate[ibody][2])*dy; + + dx = x[i][0]; + dy = x[i][1]; + dz = x[i][2]; + x[i][0] = cos(rotate[ibody][1])*dx + sin(rotate[ibody][1])*dz; + x[i][2] = -sin(rotate[ibody][1])*dx + cos(rotate[ibody][1])*dz; + + dx = x[i][0]; + dy = x[i][1]; + dz = x[i][2]; + x[i][1] = cos(rotate[ibody][0])*dy - sin(rotate[ibody][0])*dz; + x[i][2] = sin(rotate[ibody][0])*dy + cos(rotate[ibody][0])*dz; + + v[i][0] = (omega[ibody][1]*x[i][2] - omega[ibody][2]*x[i][1]) + vcm[ibody][0]; + v[i][1] = (omega[ibody][2]*x[i][0] - omega[ibody][0]*x[i][2]) + vcm[ibody][1]; + v[i][2] = (omega[ibody][0]*x[i][1] - omega[ibody][1]*x[i][0]) + vcm[ibody][2]; + + // add center of mass to displacement + // map back into periodic box via xbox,ybox,zbox + // for triclinic, add in box tilt factors as well + x[i][0] += xcm[ibody][0] - xbox*xprd; + x[i][1] += xcm[ibody][1] - ybox*yprd; + x[i][2] += xcm[ibody][2] - zbox*zprd; + + // virial = unwrapped coords dotted into body constraint force + // body constraint force = implied force due to v change minus f external + // assume f does not include forces internal to body + // 1/2 factor b/c final_integrate contributes other half + // assume per-atom contribution is due to constraint force on that atom + + if (evflag) { + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + fc0 = massone*(v[i][0] - v0)/dtf - f[i][0] + Gamma_MD[ibody]*(v0-up[i][0]); + fc1 = massone*(v[i][1] - v1)/dtf - f[i][1] + Gamma_MD[ibody]*(v1-up[i][1]); + fc2 = massone*(v[i][2] - v2)/dtf - f[i][2] + Gamma_MD[ibody]*(v2-up[i][2]); + + vr[0] = 0.5*x0*fc0; + vr[1] = 0.5*x1*fc1; + vr[2] = 0.5*x2*fc2; + vr[3] = 0.5*x0*fc1; + vr[4] = 0.5*x0*fc2; + vr[5] = 0.5*x1*fc2; + + v_tally(1,&i,1.0,vr); + } + + } + + +} +/* ---------------------------------------------------------------------- + remap xcm of each rigid body back into periodic simulation box + done during pre_neighbor so will be after call to pbc() + and after fix_deform::pre_exchange() may have flipped box + use domain->remap() in case xcm is far away from box + due to 1st definition of rigid body or due to box flip + if don't do this, then atoms of a body which drifts far away + from a triclinic box will be remapped back into box + with huge displacements when the box tilt changes via set_x() +------------------------------------------------------------------------- */ + +void FixLbRigidPCSphere::pre_neighbor() +{ + int original,oldimage,newimage; + + for (int ibody = 0; ibody < nbody; ibody++) { + original = imagebody[ibody]; + domain->remap(xcm[ibody],imagebody[ibody]); + + if (original == imagebody[ibody]) remapflag[ibody][3] = 0; + else { + oldimage = original & IMGMASK; + newimage = imagebody[ibody] & IMGMASK; + remapflag[ibody][0] = newimage - oldimage; + oldimage = (original >> IMGBITS) & IMGMASK; + newimage = (imagebody[ibody] >> IMGBITS) & IMGMASK; + remapflag[ibody][1] = newimage - oldimage; + oldimage = original >> IMG2BITS; + newimage = imagebody[ibody] >> IMG2BITS; + remapflag[ibody][2] = newimage - oldimage; + remapflag[ibody][3] = 1; + } + } + + // adjust image flags of any atom in a rigid body whose xcm was remapped + + int *atomimage = atom->image; + int nlocal = atom->nlocal; + + int ibody,idim,otherdims; + + for (int i = 0; i < nlocal; i++) { + if (body[i] == -1) continue; + if (remapflag[body[i]][3] == 0) continue; + ibody = body[i]; + + if (remapflag[ibody][0]) { + idim = atomimage[i] & IMGMASK; + otherdims = atomimage[i] ^ idim; + idim -= remapflag[ibody][0]; + idim &= IMGMASK; + atomimage[i] = otherdims | idim; + } + if (remapflag[ibody][1]) { + idim = (atomimage[i] >> IMGBITS) & IMGMASK; + otherdims = atomimage[i] ^ (idim << IMGBITS); + idim -= remapflag[ibody][1]; + idim &= IMGMASK; + atomimage[i] = otherdims | (idim << IMGBITS); + } + if (remapflag[ibody][2]) { + idim = atomimage[i] >> IMG2BITS; + otherdims = atomimage[i] ^ (idim << IMG2BITS); + idim -= remapflag[ibody][2]; + idim &= IMGMASK; + atomimage[i] = otherdims | (idim << IMG2BITS); + } + } +} +/* ---------------------------------------------------------------------- + count # of degrees-of-freedom removed by fix_rigid for atoms in igroup +------------------------------------------------------------------------- */ + +int FixLbRigidPCSphere::dof(int igroup) +{ + int groupbit = group->bitmask[igroup]; + + // ncount = # of atoms in each rigid body that are also in group + + int *mask = atom->mask; + int nlocal = atom->nlocal; + + int *ncount = new int[nbody]; + for (int ibody = 0; ibody < nbody; ibody++) ncount[ibody] = 0; + + for (int i = 0; i < nlocal; i++) + if (body[i] >= 0 && mask[i] & groupbit) ncount[body[i]]++; + + int *nall = new int[nbody]; + MPI_Allreduce(ncount,nall,nbody,MPI_INT,MPI_SUM,world); + + // remove 3N - 6 dof for each rigid body if more than 2 atoms in igroup + // remove 3N - 5 dof for each diatomic rigid body in igroup + + int n = 0; + for (int ibody = 0; ibody < nbody; ibody++) { + if (nall[ibody] > 2) n += 3*nall[ibody] - 6; + else if (nall[ibody] == 2) n++; + } + + delete [] ncount; + delete [] nall; + + return n; +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based arrays +------------------------------------------------------------------------- */ + +double FixLbRigidPCSphere::memory_usage() +{ + int nmax = atom->nmax; + double bytes = nmax * sizeof(int); + bytes += nmax*3 * sizeof(double); + bytes += maxvatom*6 * sizeof(double); + + return bytes; +} + +/* ---------------------------------------------------------------------- + allocate local atom-based arrays +------------------------------------------------------------------------- */ + +void FixLbRigidPCSphere::grow_arrays(int nmax) +{ + + memory->grow(body,nmax,"rigid:body"); + memory->grow(up,nmax,3,"rigid:up"); + +} + +/* ---------------------------------------------------------------------- + copy values within local atom-based arrays +------------------------------------------------------------------------- */ + +void FixLbRigidPCSphere::copy_arrays(int i, int j, int delflag) +{ + body[j] = body[i]; + up[j][0] = up[i][0]; + up[j][1] = up[i][1]; + up[j][2] = up[i][2]; +} + +/* ---------------------------------------------------------------------- + initialize one atom's array values, called when atom is created +------------------------------------------------------------------------- */ + +void FixLbRigidPCSphere::set_arrays(int i) +{ + body[i] = -1; + up[i][0] = 0.0; + up[i][1] = 0.0; + up[i][2] = 0.0; + +} + +/* ---------------------------------------------------------------------- + pack values in local atom-based arrays for exchange with another proc +------------------------------------------------------------------------- */ + +int FixLbRigidPCSphere::pack_exchange(int i, double *buf) +{ + buf[0] = body[i]; + buf[1] = up[i][0]; + buf[2] = up[i][1]; + buf[3] = up[i][2]; + return 4; + +} + +/* ---------------------------------------------------------------------- + unpack values in local atom-based arrays from exchange with another proc +------------------------------------------------------------------------- */ + +int FixLbRigidPCSphere::unpack_exchange(int nlocal, double *buf) +{ + body[nlocal] = static_cast (buf[0]); + up[nlocal][0] = buf[1]; + up[nlocal][1] = buf[2]; + up[nlocal][2] = buf[3]; + return 4; +} + +/* ---------------------------------------------------------------------- */ + +void FixLbRigidPCSphere::reset_dt() +{ + dtv = update->dt; + dtf = 0.5 * update->dt * force->ftm2v; +} + +/* ---------------------------------------------------------------------- + return temperature of collection of rigid bodies + non-active DOF are removed by fflag/tflag and in tfactor +------------------------------------------------------------------------- */ + +double FixLbRigidPCSphere::compute_scalar() +{ + + double inertia; + double t = 0.0; + + for (int i = 0; i < nbody; i++) { + t += masstotal[i] * (fflag[i][0]*vcm[i][0]*vcm[i][0] + + fflag[i][1]*vcm[i][1]*vcm[i][1] + \ + fflag[i][2]*vcm[i][2]*vcm[i][2]); + + // wbody = angular velocity in body frame + + inertia = 2.0*masstotal[i]*sphereradius[i]*sphereradius[i]/5.0; + + t += tflag[i][0]*inertia*omega[i][0]*omega[i][0] + + tflag[i][1]*inertia*omega[i][1]*omega[i][1] + + tflag[i][2]*inertia*omega[i][2]*omega[i][2]; + } + + t *= tfactor; + return t; +} + + +/* ---------------------------------------------------------------------- + return attributes of a rigid body + 15 values per body + xcm = 0,1,2; vcm = 3,4,5; fcm = 6,7,8; torque = 9,10,11; image = 12,13,14 +------------------------------------------------------------------------- */ + +double FixLbRigidPCSphere::compute_array(int i, int j) +{ + if (j < 3) return xcm[i][j]; + if (j < 6) return vcm[i][j-3]; + if (j < 9) return (fcm[i][j-6]+fcm_fluid[i][j-6]); + if (j < 12) return (torque[i][j-9]+torque_fluid[i][j-9]); + if (j == 12) return (imagebody[i] & IMGMASK) - IMGMAX; + if (j == 13) return (imagebody[i] >> IMGBITS & IMGMASK) - IMGMAX; + return (imagebody[i] >> IMG2BITS) - IMGMAX; +} + +/* ---------------------------------------------------------------------- */ + void FixLbRigidPCSphere::compute_up(void) + { + int *mask = atom->mask; + int nlocal = atom->nlocal; + double **x = atom->x; + int i,j,k,p,m; + 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; + int trilinear_stencil = fix_lb_fluid->trilinear_stencil; + double FfP[64]; + + + for(i=0; isublo[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; + + } + } + + } diff --git a/src/USER-LB/fix_lb_rigid_pc_sphere.h b/src/USER-LB/fix_lb_rigid_pc_sphere.h new file mode 100755 index 0000000000..612a2143a6 --- /dev/null +++ b/src/USER-LB/fix_lb_rigid_pc_sphere.h @@ -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 diff --git a/src/USER-LB/fix_lb_viscous.cpp b/src/USER-LB/fix_lb_viscous.cpp new file mode 100755 index 0000000000..03beb57f3d --- /dev/null +++ b/src/USER-LB/fix_lb_viscous.cpp @@ -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; ifixnfix; 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; jone(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); +} + diff --git a/src/USER-LB/fix_lb_viscous.h b/src/USER-LB/fix_lb_viscous.h new file mode 100755 index 0000000000..7abc86549e --- /dev/null +++ b/src/USER-LB/fix_lb_viscous.h @@ -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 diff --git a/src/special.cpp b/src/special.cpp index d455ca5cb6..170115521f 100644 --- a/src/special.cpp +++ b/src/special.cpp @@ -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"); diff --git a/src/update.cpp b/src/update.cpp index ad58842937..ab18259771 100644 --- a/src/update.cpp +++ b/src/update.cpp @@ -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;