Now we can delete particles in some region. (still have problems for …

…parallel computing)

interfaceToLammps/library.cpp

@@ -20,12 +20,15 @@
 #include "update.h"
 #include "input.h"
 #include "atom.h"
+#include "atom_vec.h"
 #include "fix_fluid_drag.h" // added JS
 #include "modify.h"
 #include "string.h"
 #include "stdio.h"
 #include "math.h"
 #include "math_const.h"
+#include "random_park.h"
+#include "memory.h"
 
 using namespace LAMMPS_NS;
 
@@ -372,3 +375,149 @@ void lammps_set_timestep(void *ptr, double dt_i)
  lammps->update->dt = dt_i;
 }
 
+/* ---------------------------------------------------------------------- */
+
+void lammps_create_particle(void *ptr)
+{
+ LAMMPS *lammps = (LAMMPS *) ptr;
+
+ int natom = static_cast<int> (lammps->atom->natoms);
+ int ntype = 1;
+ double xNew[3];
+
+ xNew[0] = 0.1;
+ xNew[1] = 0.1;
+ xNew[2] = 0.1;
+
+ lammps->atom->avec->create_atom(ntype,xNew);
+ tagint *tag = lammps->atom->tag;
+ int nlocal = lammps->atom->nlocal;
+
+ tagint max = 0;
+ tagint maxtag_all = 0;
+
+ for (int i = 0; i < nlocal; i++) max = MAX(max,tag[i]);
+ MPI_Allreduce(&max,&maxtag_all,1,MPI_LMP_TAGINT,MPI_MAX,lammps->world);
+
+
+ int n = lammps->atom->nlocal - 1;
+ lammps->atom->tag[n] = maxtag_all + 1;
+
+ lammps->atom->v[n][0] = 0.1;
+ lammps->atom->v[n][1] = 0.1;
+ lammps->atom->v[n][2] = 0.1;
+
+ double radtmp = 0.005;
+ lammps->atom->radius[n] = radtmp;
+ lammps->atom->rmass[n] = 4.0*3.14159265358917323846/3.0 * radtmp*radtmp*radtmp * 2500;
+
+ for (int j = 0; j < lammps->modify->nfix; j++)
+ if (lammps->modify->fix[j]->create_attribute) lammps->modify->fix[j]->set_arrays(n);
+
+ lammps->atom->natoms += 1;
+}
+
+
+void lammps_delete_particle(void *ptr, int* deleteList)
+{
+ LAMMPS *lammps = (LAMMPS *) ptr;
+ class RanPark *random = new RanPark(lammps,100);
+
+ int i,j,m,iwhichglobal,iwhichlocal;
+ int ndel,ndeltopo[4];
+ int *list,*mark;
+
+ // if (update->ntimestep != next_reneighbor) return;
+
+ // grow list and mark arrays if necessary
+
+ int nmax = 0;
+ if (lammps->atom->nlocal > nmax) {
+ lammps->memory->destroy(list);
+ lammps->memory->destroy(mark);
+ nmax = lammps->atom->nmax;
+ lammps->memory->create(list,nmax,"evaporate:list");
+ lammps->memory->create(mark,nmax,"evaporate:mark");
+ }
+
+ // ncount = # of deletable atoms in region that I own
+ // nall = # on all procs // nbefore = # on procs before me
+ // list[ncount] = list of local indices of atoms I can delete
+
+ double **x = lammps->atom->x;
+ int *mask = lammps->atom->mask;
+ tagint *tag = lammps->atom->tag;
+ int nlocal = lammps->atom->nlocal;
+
+ printf("listing particles..");
+ int ncount = 0;
+ for (i = 0; i < nlocal; i++)
+ {
+ printf("++++=> checking particle: %5d\n", tag[i]);
+ printf("++++=> deteleList[tag[i]-1] is: %d\n", deleteList[tag[i]-1]);
+ if (deleteList[tag[i]-1] == 1) // delete the particle when assigned in openfoam
+ list[ncount++] = i;
+ }
+
+ int nall,nbefore;
+ MPI_Allreduce(&ncount,&nall,1,MPI_INT,MPI_SUM,lammps->world);
+ MPI_Scan(&ncount,&nbefore,1,MPI_INT,MPI_SUM,lammps->world);
+ nbefore -= ncount;
+
+ // ndel = total # of atom deletions, in or out of region
+ // ndeltopo[1,2,3,4] = ditto for bonds, angles, dihedrals, impropers
+ // mark[] = 1 if deleted
+
+ ndel = 0;
+
+ for (i = 0; i < nlocal; i++) mark[i] = 0;
+
+ printf("marking particles..");
+ while (nall && ndel < 0.5) {
+ iwhichglobal = static_cast<int> (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--;
+ }
+ ndel++;
+ nall--;
+ }
+
+ // atomic deletions
+ // choose atoms randomly across all procs and mark them for deletion
+ // shrink eligible list as my atoms get marked
+ // keep ndel,ncount,nall,nbefore current after each atom deletion
+
+ // delete my marked atoms
+ // loop in reverse order to avoid copying marked atoms
+
+ printf("manupulating..");
+ AtomVec *avec = lammps->atom->avec;
+
+ for (i = nlocal-1; i >= 0; i--) {
+ if (mark[i]) {
+ avec->copy(lammps->atom->nlocal-1,i,1);
+ lammps->atom->nlocal--;
+ }
+ }
+
+ // reset global natoms and bonds, angles, etc
+ // if global map exists, reset it now instead of waiting for comm
+ // since deleting atoms messes up ghosts
+
+ lammps->atom->natoms -= ndel;
+
+ // if (ndel && atom->map_style) {
+ // atom->nghost = 0;
+ // atom->map_init();
+ // atom->map_set();
+ // }
+
+ // // statistics
+
+ // ndeleted += ndel;
+ // next_reneighbor = update->ntimestep + nevery;
+}

interfaceToLammps/library.h

@@ -55,6 +55,8 @@ extern "C" {
  void lammps_step(void* ptr, int n);
  void lammps_set_timestep(void* ptr, double dt_i);
  double lammps_get_timestep(void* ptr);
+ void lammps_create_particle(void* ptr);
+ void lammps_delete_particle(void* ptr, int* deleteList);
 
  /* used in the sorting part when assigning data from OpenFOAM*/
  struct tagpair {

lammpsFoam/softParticleCloud.C

@@ -188,7 +188,6 @@ void softParticleCloud::initLammps()
  Info << "finished moving..." << endl;
  Info<< "execution time is: " << runTime_.elapsedCpuTime() << endl;
 
-
  lammps_step(lmp_, 0);
  Info << "finished steping..." << endl;
 
@@ -340,6 +339,41 @@ void softParticleCloud::initConstructParticles
  }
 }
 
+// Add new particles in OpenFOAM
+void softParticleCloud::addNewParticles()
+{
+ Pout << "Adding new particle... " << endl;
+
+ vector pos = mesh_.C()[0];
+ label cellI = 0;
+
+ vector velo = vector(0,0,0);
+
+ scalar ds = 0.005;
+ scalar rhos = 2500;
+ label tags = 2;
+ label lmpCpuIds = 0;
+ label types = 1;
+
+ // create a new softParticle when it is in the current processor
+ // but the computer is running much slower than before.
+ softParticle* ptr =
+ new softParticle
+ (
+ pMesh(),
+ pos,
+ cellI,
+ ds,
+ velo,
+ rhos,
+ tags,
+ lmpCpuIds,
+ types
+ );
+
+ addParticle(ptr);
+}
+
 
 // Wrap up Lammps, i.e. delete the pointer.
 void softParticleCloud::finishLammps()
@@ -432,6 +466,12 @@ softParticleCloud::softParticleCloud
  cpuTimeSplit_(6, 0.0)
 {
  subCycles_ = readScalar(cloudProperties_.lookup("subCycles"));
+ addParticleFlag_ = cloudProperties_.lookupOrDefault("addParticle",0);
+ deleteParticleFlag_ = cloudProperties_.lookupOrDefault("deleteParticle",0);
+ addParticleBox_ =
+ cloudProperties_.lookupOrDefault("addParticleBox",symmTensor::zero);
+ deleteParticleBox_ =
+ cloudProperties_.lookupOrDefault("deleteParticleBox",symmTensor::zero);
 
  nLocal_ = 0;
 
@@ -507,6 +547,84 @@ void softParticleCloud::setPositionCell()
  }
 }
 
+//- Adding and deleting particles
+void softParticleCloud::addAndDeleteParticle()
+{
+ // Try adding new particles
+ if (addParticleFlag_ == 1)
+ {
+ addNewParticles();
+ lammps_create_particle(lmp_);
+ }
+
+ // Try deleting existing particles
+ if (deleteParticleFlag_ == 1)
+ {
+ int maxTag = 0;
+ int i = 0;
+ for
+ (
+ softParticleCloud::iterator pIter = begin();
+ pIter != end();
+ ++pIter, ++i
+ )
+ {
+ softParticle& p = pIter();
+ maxTag = max(p.ptag(),maxTag);
+ }
+
+ int* deleteList = new int [maxTag];
+
+ for (int i = 0; i < maxTag; i++)
+ {
+ // TODO: this may have problem for parallel computing
+ deleteList[i] = 0;
+ }
+
+ i = 0;
+ for
+ (
+ softParticleCloud::iterator pIter = begin();
+ pIter != end();
+ ++pIter, ++i
+ )
+ {
+ softParticle& p = pIter();
+
+ scalar x1 = deleteParticleBox_.component(0);
+ scalar x2 = deleteParticleBox_.component(1);
+ scalar y1 = deleteParticleBox_.component(2);
+ scalar y2 = deleteParticleBox_.component(3);
+ scalar z1 = deleteParticleBox_.component(4);
+ scalar z2 = deleteParticleBox_.component(5);
+
+ vector pPosition = p.position();
+
+ Info<< "position is: " << pPosition << endl;
+ Info<< "p tag is: " << p.ptag() << endl;
+ Info<< "deleteList is: " << deleteList[p.ptag() - 1] << endl;
+ Info<< "deleteParticleBox_ is: " << deleteParticleBox_ << endl;
+ Info<< "x1 " << (pPosition.x() - x1)*(pPosition.x() - x2) << endl;
+ Info<< "x2 " << (pPosition.y() - y1)*(pPosition.y() - y2) << endl;
+ Info<< "x3 " << (pPosition.z() - z1)*(pPosition.z() - z2) << endl;
+
+ if ((pPosition.x() - x1)*(pPosition.x() - x2) < 0 && 
+ (pPosition.y() - y1)*(pPosition.y() - y2) < 0 &&
+ (pPosition.z() - z1)*(pPosition.z() - z2) < 0 ) 
+ {
+ Info<< "deleting particle at: " << pPosition << endl; 
+ // TODO: this may have problem for parallel computing
+ deleteList[p.ptag() - 1] = 1;
+ deleteParticle(p);
+ }
+ }
+
+ lammps_delete_particle(lmp_, deleteList);
+
+ delete [] deleteList;
+ }
+}
+
 // Temp placement;
 
 // transpose information on processors
@@ -733,22 +851,33 @@ void softParticleCloud::lammpsEvolveForward
 
  cpuTimeSplit_[3] += runTime_.elapsedCpuTime() - t0;
  t0 = runTime_.elapsedCpuTime();
+
+
+ // Adding and deleting particle at certain regions
+ addAndDeleteParticle();
+
  // Ask lammps to move certain steps forward
+ Info<< "LAMMPS evolving.. " << endl;
  lammps_step(lmp_, nstep);
 
+ Info<< "finished moving the particles in LAMMPS." << endl;
  cpuTimeSplit_[4] += runTime_.elapsedCpuTime() - t0;
  t0 = runTime_.elapsedCpuTime();
  // Start getting information from LAMMPS
  // Harvest the number of particles in each lmp cpu
  int lmpNLocal = lammps_get_local_n(lmp_);
 
+ Info<< "the number of particles in LAMMPS now is: " << lmpNLocal << endl;
+
  // Harvest more infomation from each lmp cpu
  double* fromLmpXArrayLocal = new double [3*lmpNLocal];
  double* fromLmpVArrayLocal = new double [3*lmpNLocal];
  int* fromLmpFoamCpuIdArrayLocal = new int[lmpNLocal];
  int* fromLmpLmpCpuIdArrayLocal = new int[lmpNLocal];
  int* fromLmpTagArrayLocal = new int[lmpNLocal];
 
+ Info<< "getting local info.." << endl;
+
  lammps_get_local_info
  (
  lmp_,
@@ -759,6 +888,8 @@ void softParticleCloud::lammpsEvolveForward
  fromLmpTagArrayLocal
  );
 
+ Info<< "local info obtained!" << endl;
+
  // Transform the data obtained from lammps to openfoam format
  vectorList fromLmpXList(lmpNLocal, vector::zero);
  vectorList fromLmpVList(lmpNLocal, vector::zero);
@@ -786,12 +917,20 @@ void softParticleCloud::lammpsEvolveForward
  fromLmpVArrayLocal[3*i + 2]
  );
 
+ Info<< "fromLmpXArrayLocal is: " << fromLmpXArrayLocal[i] << endl;
+ Info<< "fromLmpVArrayLocal is: " << fromLmpVArrayLocal[i] << endl;
+ Info<< "fromLmpFoamCpuIdArrayLocal is: " << fromLmpFoamCpuIdArrayLocal[i] << endl;
+ Info<< "fromLmpLmpCpuIdArrayLocal is: " << fromLmpLmpCpuIdArrayLocal[i] << endl;
+ Info<< "fromLmpTagArrayLocal is: " << fromLmpTagArrayLocal[i] << endl;
+
  label foamCpuId = fromLmpFoamCpuIdArrayLocal[i];
  fromLmpFoamCpuIdList[i] = foamCpuId;
  fromLmpTagList[i] = fromLmpTagArrayLocal[i];
  foamParticleNo[foamCpuId] ++;
  }
 
+ Info<< "data transformed into each processor!" << endl;
+
  delete [] fromLmpXArrayLocal;
  delete [] fromLmpVArrayLocal;
  delete [] fromLmpFoamCpuIdArrayLocal;