From d32f0a53f2fc0ea9da92282204f592cce6473758 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 18 Mar 2009 17:37:46 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@2676 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/compute_displace_atom.cpp | 2 +- src/compute_pe_atom.cpp | 1 + src/compute_stress_atom.cpp | 1 + src/fix_evaporate.cpp | 191 ++++++++++++++++++++++++++++++++++ src/fix_evaporate.h | 42 ++++++++ src/read_data.cpp | 2 +- src/style.h | 2 + 7 files changed, 239 insertions(+), 2 deletions(-) create mode 100644 src/fix_evaporate.cpp create mode 100644 src/fix_evaporate.h diff --git a/src/compute_displace_atom.cpp b/src/compute_displace_atom.cpp index 62b41cbfa8..5174e80843 100644 --- a/src/compute_displace_atom.cpp +++ b/src/compute_displace_atom.cpp @@ -78,7 +78,7 @@ void ComputeDisplaceAtom::compute_peratom() // grow local displacement array if necessary - if (atom->nmax > nmax) { + if (atom->nlocal > nmax) { memory->destroy_2d_double_array(displace); nmax = atom->nmax; displace = diff --git a/src/compute_pe_atom.cpp b/src/compute_pe_atom.cpp index 31ac624480..87042ec430 100755 --- a/src/compute_pe_atom.cpp +++ b/src/compute_pe_atom.cpp @@ -80,6 +80,7 @@ void ComputePEAtom::compute_peratom() error->all("Per-atom energy was not tallied on needed timestep"); // grow local energy array if necessary + // needs to be atom->nmax in length if (atom->nmax > nmax) { memory->sfree(energy); diff --git a/src/compute_stress_atom.cpp b/src/compute_stress_atom.cpp index 08bfb03164..cdb13f3116 100644 --- a/src/compute_stress_atom.cpp +++ b/src/compute_stress_atom.cpp @@ -94,6 +94,7 @@ void ComputeStressAtom::compute_peratom() error->all("Per-atom virial was not tallied on needed timestep"); // grow local stress array if necessary + // needs to be atom->nmax in length if (atom->nmax > nmax) { memory->destroy_2d_double_array(stress); diff --git a/src/fix_evaporate.cpp b/src/fix_evaporate.cpp new file mode 100644 index 0000000000..f74d18c9ad --- /dev/null +++ b/src/fix_evaporate.cpp @@ -0,0 +1,191 @@ +/* ---------------------------------------------------------------------- + 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 "stdlib.h" +#include "string.h" +#include "fix_evaporate.h" +#include "atom.h" +#include "atom_vec.h" +#include "update.h" +#include "domain.h" +#include "region.h" +#include "comm.h" +#include "random_park.h" +#include "random_mars.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define MIN(A,B) ((A) < (B)) ? (A) : (B) +#define MAX(A,B) ((A) > (B)) ? (A) : (B) + +/* ---------------------------------------------------------------------- */ + +FixEvaporate::FixEvaporate(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + if (narg != 7) error->all("Illegal fix evaporate command"); + + scalar_flag = 1; + scalar_vector_freq = 1; + extscalar = 0; + + nevery = atoi(arg[3]); + nflux = atoi(arg[4]); + iregion = domain->find_region(arg[5]); + int seed = atoi(arg[6]); + + if (nevery <= 0 || nflux <= 0) error->all("Illegal fix evaporate command"); + if (iregion == -1) error->all("Fix evaporate region ID does not exist"); + if (seed <= 0) error->all("Illegal fix evaporate command"); + + // random number generator, same for all procs + + random = new RanPark(lmp,seed); + + // set up reneighboring + + force_reneighbor = 1; + next_reneighbor = (update->ntimestep/nevery)*nevery + nevery; + ndeleted = 0; + + nmax = 0; + list = NULL; + mark = NULL; +} + +/* ---------------------------------------------------------------------- */ + +FixEvaporate::~FixEvaporate() +{ + delete random; + memory->sfree(list); + memory->sfree(mark); +} + +/* ---------------------------------------------------------------------- */ + +int FixEvaporate::setmask() +{ + int mask = 0; + mask |= PRE_EXCHANGE; + return mask; +} + +/* ---------------------------------------------------------------------- + perform particle deletion + done before exchange, borders, reneighbor + so that ghost atoms and neighbor lists will be correct +------------------------------------------------------------------------- */ + +void FixEvaporate::pre_exchange() +{ + int i,iwhichglobal,iwhichlocal; + + if (update->ntimestep != next_reneighbor) return; + + // grow list and mark arrays if necessary + + if (atom->nlocal > nmax) { + memory->sfree(list); + memory->sfree(mark); + nmax = atom->nmax; + list = (int *) memory->smalloc(nmax*sizeof(int),"fix/evaporate:list"); + mark = (int *) memory->smalloc(nmax*sizeof(int),"fix/evaporate:mark"); + } + + // nall = # of deletable atoms in region + // nbefore = # on procs before me + + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + int ncount = 0; + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + if (domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2])) + list[ncount++] = i; + + int nall,nbefore; + MPI_Allreduce(&ncount,&nall,1,MPI_INT,MPI_SUM,world); + MPI_Scan(&ncount,&nbefore,1,MPI_INT,MPI_SUM,world); + nbefore -= ncount; + + // nwhack = total number of atoms to delete + // choose atoms randomly across all procs and mark them for deletion + // shrink local list of candidates as my atoms get marked for deletion + + int nwhack = MIN(nflux,nall); + + for (i = 0; i < nlocal; i++) mark[i] = 0; + + for (i = 0; i < nwhack; i++) { + iwhichglobal = static_cast (nall*random->uniform()); + if (iwhichglobal < nbefore) nbefore--; + else if (iwhichglobal < nbefore + ncount) { + iwhichlocal = iwhichglobal - nbefore; + mark[list[iwhichlocal]] = 1; + list[iwhichlocal] = list[ncount-1]; + ncount--; + } + nall--; + } + + // delete my marked atoms + // loop in reverse order to avoid copying marked atoms + + AtomVec *avec = atom->avec; + + for (i = nlocal-1; i >= 0; i--) { + if (mark[i]) { + avec->copy(atom->nlocal-1,i); + atom->nlocal--; + } + } + + // reset global natoms + // if global map exists, reset it + + atom->natoms -= nwhack; + if (nwhack && atom->map_style) { + atom->map_init(); + atom->map_set(); + } + + // statistics + + ndeleted += nwhack; + next_reneighbor = update->ntimestep + nevery; +} + +/* ---------------------------------------------------------------------- + return number of deleted particles +------------------------------------------------------------------------- */ + +double FixEvaporate::compute_scalar() +{ + return 1.0*ndeleted; +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based arrays +------------------------------------------------------------------------- */ + +double FixEvaporate::memory_usage() +{ + double bytes = 2*nmax * sizeof(int); + return bytes; +} diff --git a/src/fix_evaporate.h b/src/fix_evaporate.h new file mode 100644 index 0000000000..5b9866222e --- /dev/null +++ b/src/fix_evaporate.h @@ -0,0 +1,42 @@ +/* ---------------------------------------------------------------------- + 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 FIX_EVAPORATE_H +#define FIX_EVAPORATE_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixEvaporate : public Fix { + public: + FixEvaporate(class LAMMPS *, int, char **); + ~FixEvaporate(); + int setmask(); + void pre_exchange(); + double compute_scalar(); + double memory_usage(); + + private: + int nevery,nflux,iregion; + int ndeleted; + + int nmax; + int *list,*mark; + + class RanPark *random; +}; + +} + +#endif diff --git a/src/read_data.cpp b/src/read_data.cpp index a2ae8cef8c..9673e94486 100644 --- a/src/read_data.cpp +++ b/src/read_data.cpp @@ -441,7 +441,7 @@ void ReadData::atoms() if (natoms != atom->natoms) error->all("Did not assign all atoms correctly"); // if any atom ID < 0, error - // if all atom IDs == 0, tag_enable = 0 + // if all atom IDs = 0, tag_enable = 0 // if any atom ID > 0, error if any atom ID == 0 // not checking if atom IDs > natoms or are unique diff --git a/src/style.h b/src/style.h index ada811be8b..c6bc3fecbb 100644 --- a/src/style.h +++ b/src/style.h @@ -156,6 +156,7 @@ DumpStyle(xyz,DumpXYZ) #include "fix_dt_reset.h" #include "fix_efield.h" #include "fix_enforce2d.h" +#include "fix_evaporate.h" #include "fix_gravity.h" #include "fix_gyration.h" #include "fix_heat.h" @@ -216,6 +217,7 @@ FixStyle(drag,FixDrag) FixStyle(dt/reset,FixDtReset) FixStyle(efield,FixEfield) FixStyle(enforce2d,FixEnforce2D) +FixStyle(evaporate,FixEvaporate) FixStyle(gravity,FixGravity) FixStyle(gyration,FixGyration) FixStyle(heat,FixHeat)