From 61c291ba7b90ab87adef6bcfc7f04dfe6539c6a8 Mon Sep 17 00:00:00 2001 From: iurnus Date: Fri, 23 Jan 2015 09:59:13 -0500 Subject: [PATCH] Now we can delete particles in some region. (still have problems for parallel computing) --- interfaceToLammps/library.cpp | 149 +++++++++++++++++++++++++++++++++ interfaceToLammps/library.h | 2 + lammpsFoam/softParticleCloud.C | 141 ++++++++++++++++++++++++++++++- lammpsFoam/softParticleCloud.H | 14 ++++ 4 files changed, 305 insertions(+), 1 deletion(-) diff --git a/interfaceToLammps/library.cpp b/interfaceToLammps/library.cpp index 38abbe9..28ce416 100644 --- a/interfaceToLammps/library.cpp +++ b/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 (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 (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; +} diff --git a/interfaceToLammps/library.h b/interfaceToLammps/library.h index ef53cfa..f4adead 100644 --- a/interfaceToLammps/library.h +++ b/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 { diff --git a/lammpsFoam/softParticleCloud.C b/lammpsFoam/softParticleCloud.C index f10b861..f8b9402 100644 --- a/lammpsFoam/softParticleCloud.C +++ b/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,15 +851,24 @@ 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]; @@ -749,6 +876,8 @@ void softParticleCloud::lammpsEvolveForward 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; diff --git a/lammpsFoam/softParticleCloud.H b/lammpsFoam/softParticleCloud.H index 7d2bad1..da9aad2 100644 --- a/lammpsFoam/softParticleCloud.H +++ b/lammpsFoam/softParticleCloud.H @@ -156,6 +156,9 @@ class softParticleCloud int* type ); + //- Add OpenFOAM particles + void addNewParticles(); + //- Disallow default bitwise copy construct softParticleCloud(const softParticleCloud&); @@ -187,6 +190,15 @@ protected: IOdictionary cloudProperties_; + //- If add any particles in the simulation + scalar addParticleFlag_; + symmTensor addParticleBox_; + + + //- If delete any particles in the simulation + scalar deleteParticleFlag_; + symmTensor deleteParticleBox_; + // References to the physical field //- Velocity @@ -224,6 +236,8 @@ protected: // move across the processor boundary void setPositionCell(); + //- Add and delete particles + void addAndDeleteParticle(); //- Transpose the scalar of different processors template