diff --git a/src/DSMC/Install.csh b/src/DSMC/Install.csh new file mode 100644 index 0000000000..2b0f0fdc38 --- /dev/null +++ b/src/DSMC/Install.csh @@ -0,0 +1,20 @@ +# Install/unInstall package classes in LAMMPS + +if ($1 == 1) then + + cp style_dsmc.h .. + + cp pair_dsmc.cpp .. + + cp pair_dsmc.h .. + +else if ($1 == 0) then + + rm ../style_dsmc.h + touch ../style_dsmc.h + + rm ../pair_dsmc.cpp + + rm ../pair_dsmc.h + +endif diff --git a/src/DSMC/pair_dsmc.cpp b/src/DSMC/pair_dsmc.cpp new file mode 100644 index 0000000000..dc97f7fb85 --- /dev/null +++ b/src/DSMC/pair_dsmc.cpp @@ -0,0 +1,525 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_dsmc.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "memory.h" +#include "error.h" +#include "domain.h" +#include "update.h" +#include "random_mars.h" +#include "limits.h" + +using namespace LAMMPS_NS; + +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#define MAX(a,b) ((a) > (b) ? (a) : (b)) + +/* ---------------------------------------------------------------------- */ + +PairDSMC::PairDSMC(LAMMPS *lmp) : Pair(lmp) +{ + single_enable = 0; + + total_number_of_collisions = 0; + max_particles = max_particle_list = 0; + next_particle = NULL; + random = NULL; +} + +/* ---------------------------------------------------------------------- */ + +PairDSMC::~PairDSMC() +{ + if (allocated) { + memory->destroy_2d_int_array(setflag); + memory->destroy_2d_double_array(cutsq); + memory->destroy_2d_double_array(sigma); + memory->destroy_2d_double_array(cut); + memory->destroy_2d_double_array(V_sigma_max); + memory->destroy_2d_int_array(particle_list); + memory->destroy_2d_int_array(first); + memory->destroy_2d_int_array(number); + } + + delete [] next_particle; + delete random; +} + +/* ---------------------------------------------------------------------- */ + +void PairDSMC::compute(int eflag, int vflag) +{ + double **x = atom->x; + double *mass = atom->mass; + int *type = atom->type; + int nlocal = atom->nlocal; + + for (int i = 1; i <= atom->ntypes; ++i) + for (int j = 0; j < total_ncells; ++j) { + first[i][j] = -1; + number[i][j] = 0; + } + + if (atom->nmax > max_particles) { + delete [] next_particle; + max_particles = atom->nmax; + next_particle = new int[max_particles]; + } + + // find each particle's cell and sort by type + // assume a constant volume and shape simulation domain + // skip particle if outside processor domain + + for (int i = 0; i < nlocal; ++i) { + int xcell = static_cast((x[i][0] - domain->boxlo[0])/cellx); + int ycell = static_cast((x[i][1] - domain->boxlo[1])/celly); + int zcell = static_cast((x[i][2] - domain->boxlo[2])/cellz); + + if ((xcell < 0) or (xcell > ncellsx-1) or + (ycell < 0) or (ycell > ncellsy-1) or + (zcell < 0) or (zcell > ncellsz-1)) continue; + + int icell = xcell + ycell*ncellsx + zcell*ncellsx*ncellsy; + itype = type[i]; + next_particle[i] = first[itype][icell]; + first[itype][icell] = i; + number[itype][icell]++; + } + + for (int icell = 0; icell < total_ncells; ++icell) { + + for (itype = 1; itype <= atom->ntypes; ++itype) { + number_of_A = number[itype][icell]; + if (number_of_A > max_particle_list) { + max_particle_list = number_of_A; + particle_list = memory->grow_2d_int_array(particle_list,atom->ntypes+1, + max_particle_list, + "pair:particle_list"); + } + + int m = first[itype][icell]; + for (int k = 0; k < number_of_A; k++) { + particle_list[itype][k] = m; + m = next_particle[m]; + } + } + + for (itype = 1; itype <= atom->ntypes; ++itype) { + imass = mass[itype]; + number_of_A = number[itype][icell]; + + for (jtype = itype; jtype <= atom->ntypes; ++jtype) { + jmass = mass[jtype]; + number_of_B = number[jtype][icell]; + + reduced_mass = imass*jmass/(imass + jmass); + total_mass = imass + jmass; + jmass_tmass = jmass/total_mass; + imass_tmass = imass/total_mass; + + // if necessary, recompute V_sigma_max values + + if (recompute_vsigmamax_stride && + (update->ntimestep % recompute_vsigmamax_stride == 0)) + recompute_V_sigma_max(icell); + + // # of collisions to perform for itype-jtype pairs + + double &Vs_max = V_sigma_max[itype][jtype]; + double num_of_collisions_double = 0.5 * number_of_A * number_of_B * + weighting * Vs_max * update->dt / vol; + + if ((itype == jtype) and number_of_B) + num_of_collisions_double *= + double(number_of_B - 1) / double(number_of_B); + + int num_of_collisions = + convert_double_to_equivalent_int(num_of_collisions_double); + + if (num_of_collisions > number_of_A) + error->warning("num_of_collisions > number_of_A"); + if (num_of_collisions > number_of_B) + error->warning("num_of_collisions > number_of_B"); + + // perform collisions on pairs of particles in icell + + for (int k = 0; k < num_of_collisions; k++) { + if ((number_of_A < 1) or (number_of_B < 1)) break; + int ith_A = static_cast(random->uniform()*number_of_A); + int jth_B = static_cast(random->uniform()*number_of_B); + int i = particle_list[itype][ith_A]; + int j = particle_list[jtype][jth_B]; + if (i == j) { + k--; + continue; + } + double probability = V_sigma(i,j)/Vs_max; + if (probability > random->uniform()) scatter_random(i,j,icell); + } + } + } + } +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairDSMC::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq"); + + cut = memory->create_2d_double_array(n+1,n+1,"pair:cut"); + sigma = memory->create_2d_double_array(n+1,n+1,"pair:sigma"); + V_sigma_max = memory->create_2d_double_array(n+1,n+1,"pair:V_sigma_max"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairDSMC::settings(int narg, char **arg) +{ + if (narg != 6) error->all("Illegal pair_style command"); + + cut_global = 0.0; + max_cell_size = force->numeric(arg[0]); + seed = force->inumeric(arg[1]); + weighting = force->numeric(arg[2]); + T_ref = force->numeric(arg[3]); + recompute_vsigmamax_stride = force->inumeric(arg[4]); + vsigmamax_samples = force->inumeric(arg[5]); + + // initialize Marsaglia RNG with processor-unique seed + + if (max_cell_size <= 0.0) error->all("Illegal pair_style command"); + if (seed <= 0) error->all("Illegal pair_style command"); + if (random) delete random; + random = new RanMars(lmp,seed + comm->me); + + kT_ref = force->boltz*T_ref; + + // reset cutoffs that have been explicitly set + + if (allocated) { + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i+1; j <= atom->ntypes; j++) + if (setflag[i][j]) cut[i][j] = cut_global; + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairDSMC::coeff(int narg, char **arg) +{ + if (narg < 3 || narg > 4) error->all("Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + force->bounds(arg[0],atom->ntypes,ilo,ihi); + force->bounds(arg[1],atom->ntypes,jlo,jhi); + + double sigma_one = force->numeric(arg[2]); + + double cut_one = cut_global; + if (narg == 4) cut_one = force->numeric(arg[3]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + sigma[i][j] = sigma_one; + cut[i][j] = cut_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all("Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairDSMC::init_style() +{ + ncellsx = ncellsy = ncellsz = 1; + while (((domain->boxhi[0] - domain->boxlo[0])/ncellsx) > max_cell_size) + ncellsx++; + while (((domain->boxhi[1] - domain->boxlo[1])/ncellsy) > max_cell_size) + ncellsy++; + while (((domain->boxhi[2] - domain->boxlo[2])/ncellsz) > max_cell_size) + ncellsz++; + + cellx = (domain->boxhi[0] - domain->boxlo[0])/ncellsx; + celly = (domain->boxhi[1] - domain->boxlo[1])/ncellsy; + cellz = (domain->boxhi[2] - domain->boxlo[2])/ncellsz; + + if (comm->me == 0) { + if (screen) fprintf(screen,"DSMC cell size = %g x %g x %g\n", + cellx,celly,cellz); + if (logfile) fprintf(logfile,"DSMC cell size = %g x %g x %g\n", + cellx,celly,cellz); + } + + total_ncells = ncellsx*ncellsy*ncellsz; + vol = cellx*celly*cellz; + + particle_list = memory->create_2d_int_array(atom->ntypes+1,0, + "pair:particle_list"); + first = memory->create_2d_int_array(atom->ntypes+1,total_ncells, + "pair:first"); + number = memory->create_2d_int_array(atom->ntypes+1,total_ncells, + "pair:number"); + + for (int i = 1; i <= atom->ntypes; i++) + for (int j = 1; j <= atom->ntypes; j++) + V_sigma_max[i][j] = 0.0; + + two_pi = 8.0*atan(1.0); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairDSMC::init_one(int i, int j) +{ + if (setflag[i][j] == 0) cut[i][j] = 0.0; + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairDSMC::write_restart(FILE *fp) +{ + write_restart_settings(fp); + + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + fwrite(&setflag[i][j],sizeof(int),1,fp); + if (setflag[i][j]) { + fwrite(&sigma[i][j],sizeof(double),1,fp); + fwrite(&cut[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairDSMC::read_restart(FILE *fp) +{ + read_restart_settings(fp); + allocate(); + + int i,j; + int me = comm->me; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); + MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (setflag[i][j]) { + if (me == 0) { + fread(&sigma[i][j],sizeof(double),1,fp); + fread(&cut[i][j],sizeof(double),1,fp); + } + MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairDSMC::write_restart_settings(FILE *fp) +{ + fwrite(&cut_global,sizeof(double),1,fp); + fwrite(&max_cell_size,sizeof(double),1,fp); + fwrite(&seed,sizeof(int),1,fp); + fwrite(&offset_flag,sizeof(int),1,fp); + fwrite(&mix_flag,sizeof(int),1,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairDSMC::read_restart_settings(FILE *fp) +{ + if (comm->me == 0) { + fread(&cut_global,sizeof(double),1,fp); + fread(&max_cell_size,sizeof(double),1,fp); + fread(&seed,sizeof(int),1,fp); + fread(&offset_flag,sizeof(int),1,fp); + fread(&mix_flag,sizeof(int),1,fp); + } + + MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&max_cell_size,1,MPI_DOUBLE,0,world); + MPI_Bcast(&seed,1,MPI_INT,0,world); + MPI_Bcast(&offset_flag,1,MPI_INT,0,world); + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); + + // initialize Marsaglia RNG with processor-unique seed + // same seed that pair_style command initially specified + + if (random) delete random; + random = new RanMars(lmp,seed + comm->me); +} + +/*------------------------------------------------------------------------- + rezero and recompute the V_sigma_max values this timestep for use during + the next nrezero timesteps +-------------------------------------------------------------------------*/ + +void PairDSMC::recompute_V_sigma_max(int icell) +{ + int i,j,k; + double Vsigma_max = 0; + + if (number_of_A && number_of_B) { + for (k = 0; k < vsigmamax_samples; k++) { + i = particle_list[itype] + [static_cast(random->uniform()*number_of_A)]; + j = particle_list[jtype] + [static_cast(random->uniform()*number_of_B)]; + if (i == j) continue; + Vsigma_max = MAX(Vsigma_max,V_sigma(i,j)); + } + } + V_sigma_max[itype][jtype] = Vsigma_max; +} + +/*------------------------------------------------------------------------- + VHS model + compute the velocity vector difference between i and j and multiply by + their combined collision cross section, sigma, for neutral-neutral + collisions using the Variable Hard Sphere model +-------------------------------------------------------------------------*/ + +double PairDSMC::V_sigma(int i, int j) +{ + double relative_velocity_sq,relative_velocity,pair_sigma; + double delv[3]; + double *vi = atom->v[i]; + double *vj = atom->v[j]; + + subtract3d(vi,vj,delv); + relative_velocity_sq = dot3d(delv,delv); + relative_velocity = sqrt(relative_velocity_sq); + + // from Bird eq 4.63, and omega=0.67 + // (omega - 0.5) = 0.17 + // 1/GAMMA(2.5 - omega) = 1.06418029298371 + + if (relative_velocity_sq != 0.0) + pair_sigma = sigma[itype][jtype]* + pow(kT_ref/(0.5*reduced_mass*relative_velocity_sq),0.17) * + 1.06418029298371; + else + pair_sigma = 0.0; + + return relative_velocity*pair_sigma; +} + +/*------------------------------------------------------------------------- + generate new velocities for collided particles +-------------------------------------------------------------------------*/ + +void PairDSMC::scatter_random(int i, int j, int icell) +{ + double mag_delv,cos_phi,cos_squared,r,theta; + double delv[3],vcm[3]; + double *vi = atom->v[i]; + double *vj = atom->v[j]; + + subtract3d(vi,vj,delv); + if (itype == jtype) mag_delv = sqrt(dot3d(delv,delv))*0.5; + else mag_delv = sqrt(dot3d(delv,delv)); + + cos_phi = 1.0 - (2.0*random->uniform()); + cos_squared = MIN(1.0,cos_phi*cos_phi); + r = sqrt(1.0 - cos_squared); + delv[0] = cos_phi*mag_delv; + theta = two_pi*random->uniform(); + delv[1] = r*mag_delv*cos(theta); + delv[2] = r*mag_delv*sin(theta); + + if (itype == jtype) { + vcm[0] = (vi[0]+vj[0])*0.5; + vcm[1] = (vi[1]+vj[1])*0.5; + vcm[2] = (vi[2]+vj[2])*0.5; + vi[0] = vcm[0] + delv[0]; + vi[1] = vcm[1] + delv[1]; + vi[2] = vcm[2] + delv[2]; + vj[0] = vcm[0] - delv[0]; + vj[1] = vcm[1] - delv[1]; + vj[2] = vcm[2] - delv[2]; + } else { + vcm[0] = vi[0]*imass_tmass + vj[0]*jmass_tmass; + vcm[1] = vi[1]*imass_tmass + vj[1]*jmass_tmass; + vcm[2] = vi[2]*imass_tmass + vj[2]*jmass_tmass; + vi[0] = vcm[0] + delv[0]*jmass_tmass; + vi[1] = vcm[1] + delv[1]*jmass_tmass; + vi[2] = vcm[2] + delv[2]*jmass_tmass; + vj[0] = vcm[0] - delv[0]*imass_tmass; + vj[1] = vcm[1] - delv[1]*imass_tmass; + vj[2] = vcm[2] - delv[2]*imass_tmass; + } + + total_number_of_collisions++; +} + +/* ---------------------------------------------------------------------- + This method converts the double supplied by the calling function into + an int, which is returned. By adding a random number between 0 and 1 + to the double before converting it to an int, we ensure that, + statistically, we round down with probability identical to the + remainder and up the rest of the time. So even though we're using an + integer, we're statistically matching the exact expression represented + by the double. +------------------------------------------------------------------------- */ + +int PairDSMC::convert_double_to_equivalent_int(double input_double) +{ + if (input_double > INT_MAX) + error->all("Tried to convert a double to int, but input_double > INT_MAX"); + + int output_int = static_cast(input_double + random->uniform()); + return output_int; +} diff --git a/src/DSMC/pair_dsmc.h b/src/DSMC/pair_dsmc.h new file mode 100644 index 0000000000..eb3367d4ae --- /dev/null +++ b/src/DSMC/pair_dsmc.h @@ -0,0 +1,102 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifndef PAIR_DSMC_H +#define PAIR_DSMC_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairDSMC : public Pair { + public: + PairDSMC(class LAMMPS *); + virtual ~PairDSMC(); + virtual void compute(int, int); + virtual void settings(int, char **); + void coeff(int, char **); + void init_style(); + double init_one(int, int); + void write_restart(FILE *); + void read_restart(FILE *); + virtual void write_restart_settings(FILE *); + virtual void read_restart_settings(FILE *); + + private: + double cut_global; + double **cut; + double **sigma; + + double cellx; + double celly; + double cellz; + int ncellsx; + int ncellsy; + int ncellsz; + int total_ncells; + int total_number_of_collisions; + int recompute_vsigmamax_stride; + int vsigmamax_samples; + double T_ref; + double kT_ref; + double two_pi; + double max_cell_size; + + int seed; + int number_of_A; + int number_of_B; + int max_particle_list; + + class RanMars *random; + + int **particle_list; + int **first; + int **number; + + double **V_sigma_max; + + int max_particles; + int *next_particle; + + int itype; + int jtype; + + double imass; + double jmass; + double total_mass; + double reduced_mass; + double imass_tmass; + double jmass_tmass; + double vol; + double weighting; + + void allocate(); + void recompute_V_sigma_max(int); + double V_sigma(int, int); + void scatter_random(int, int, int); + int convert_double_to_equivalent_int(double); + + inline void subtract3d(const double *v1, const double *v2, double *v3) { + v3[0] = v2[0] - v1[0]; + v3[1] = v2[1] - v1[1]; + v3[2] = v2[2] - v1[2]; + } + + inline double dot3d(const double *v1, const double *v2) { + return( v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2] ); + } +}; + +} + +#endif diff --git a/src/DSMC/style_dsmc.h b/src/DSMC/style_dsmc.h new file mode 100644 index 0000000000..4af022dacd --- /dev/null +++ b/src/DSMC/style_dsmc.h @@ -0,0 +1,20 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PairInclude +#include "pair_dsmc.h" +#endif + +#ifdef PairClass +PairStyle(dsmc,PairDSMC) +#endif diff --git a/src/Makefile b/src/Makefile index 6f37610abb..407d571d32 100755 --- a/src/Makefile +++ b/src/Makefile @@ -13,7 +13,7 @@ OBJ = $(SRC:.cpp=.o) # Package variables -PACKAGE = asphere class2 colloid dipole dpd gpu granular \ +PACKAGE = asphere class2 colloid dipole dpd dsmc gpu granular \ kspace manybody meam molecule opt peri poems prd reax xtc PACKUSER = user-ackland user-atc user-cd-eam user-cg-cmm user-ewaldn \ diff --git a/src/style.h b/src/style.h index d35ddd3cad..9a4301a44e 100644 --- a/src/style.h +++ b/src/style.h @@ -369,26 +369,6 @@ RegionStyle(sphere,RegSphere) RegionStyle(union,RegUnion) #endif -// style files for standard packages - -#include "style_asphere.h" -#include "style_class2.h" -#include "style_colloid.h" -#include "style_dipole.h" -#include "style_dpd.h" -#include "style_gpu.h" -#include "style_granular.h" -#include "style_kspace.h" -#include "style_manybody.h" -#include "style_meam.h" -#include "style_molecule.h" -#include "style_opt.h" -#include "style_peri.h" -#include "style_poems.h" -#include "style_prd.h" -#include "style_reax.h" -#include "style_xtc.h" - // package and user add-ons #include "style_packages.h" diff --git a/src/style_packages.h b/src/style_packages.h index 1514e956cf..6120b936d4 100644 --- a/src/style_packages.h +++ b/src/style_packages.h @@ -11,13 +11,14 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -// style flies for packages +// style files for standard packages #include "style_asphere.h" #include "style_class2.h" #include "style_colloid.h" #include "style_dipole.h" #include "style_dpd.h" +#include "style_dsmc.h" #include "style_gpu.h" #include "style_granular.h" #include "style_kspace.h" @@ -27,5 +28,6 @@ #include "style_opt.h" #include "style_peri.h" #include "style_poems.h" +#include "style_prd.h" #include "style_reax.h" #include "style_xtc.h"