diff --git a/src/XTC/dump_xtc.cpp b/src/XTC/dump_xtc.cpp index 04c3fd8930..6a9cbda8e5 100644 --- a/src/XTC/dump_xtc.cpp +++ b/src/XTC/dump_xtc.cpp @@ -15,10 +15,11 @@ Contributing authors: Naveen Michaud-Agrawal (Johns Hopkins U) open-source XDR routines from Frans van Hoesel (http://md.chem.rug.nl/hoesel) - are also included in this file - Axel Kohlmeyer (UPenn) + are included in this file + Axel Kohlmeyer (Temple U) port to platforms without XDR support added support for unwrapped trajectories + support for groups ------------------------------------------------------------------------- */ #include "math.h" @@ -31,6 +32,7 @@ #include "atom.h" #include "update.h" #include "group.h" +#include "output.h" #include "error.h" #include "memory.h" @@ -52,11 +54,12 @@ int xdr3dfcoord(XDR *, float *, int *, float *); DumpXTC::DumpXTC(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg) { if (narg != 5) error->all("Illegal dump xtc command"); - if (igroup != group->find("all")) error->all("Dump xtc must use group all"); if (binary || compressed || multifile || multiproc) error->all("Invalid dump xtc filename"); - size_one = 4; + size_one = 3; + sort_flag = 1; + sortcol = 0; format_default = NULL; flush_flag = 0; unwrap_flag = 0; @@ -64,10 +67,10 @@ DumpXTC::DumpXTC(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg) // allocate global array for atom coords - natoms = static_cast (atom->natoms); + if (igroup == 0) natoms = static_cast (atom->natoms); + else natoms = static_cast (group->count(igroup)); if (natoms <= 0) error->all("Invalid natoms for dump xtc"); - if (atom->tag_consecutive() == 0) - error->all("Atom IDs must be consecutive for dump xtc"); + coords = (float *) memory->smalloc(3*natoms*sizeof(float),"dump:coords"); // sfactor = conversion of coords to XTC units @@ -77,6 +80,7 @@ DumpXTC::DumpXTC(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg) if (strcmp(update->unit_style,"lj") == 0) sfactor = 1.0; openfile(); + nevery_save = 0; ntotal = 0; } @@ -94,45 +98,26 @@ DumpXTC::~DumpXTC() /* ---------------------------------------------------------------------- */ -void DumpXTC::init() +void DumpXTC::init_style() { + if (sort_flag == 0 || sortcol != 0) + error->all("Dump xtc requires sorting by atom ID"); + // check that flush_flag is not set since dump::write() will use it if (flush_flag) error->all("Cannot set dump_modify flush for dump xtc"); -} -/* ---------------------------------------------------------------------- */ + // check that dump frequency has not changed and is not a variable -int DumpXTC::modify_param(int narg, char **arg) -{ - if (strcmp(arg[0],"unwrap") == 0) { - if (narg < 2) error->all("Illegal dump_modify command"); - if (strcmp(arg[1],"yes") == 0) unwrap_flag = 1; - else if (strcmp(arg[1],"no") == 0) unwrap_flag = 0; - else error->all("Illegal dump_modify command"); - return 2; - } else if (strcmp(arg[0],"precision") == 0) { - if (narg < 2) error->all("Illegal dump_modify command"); - precision = atof(arg[1]); - if ((fabs(precision-10.0) > EPS) && (fabs(precision-100.0) > EPS) && - (fabs(precision-1000.0) > EPS) && (fabs(precision-10000.0) > EPS) && - (fabs(precision-100000.0) > EPS) && - (fabs(precision-1000000.0) > EPS)) - error->all("Illegal dump_modify command"); - return 2; - } - return 0; -} + int idump; + for (idump = 0; idump < output->ndump; idump++) + if (strcmp(id,output->dump[idump]->id) == 0) break; + if (output->every_dump[idump] == 0) + error->all("Cannot use variable every setting for dump dcd"); -/* ---------------------------------------------------------------------- - return # of bytes of allocated memory in buf and global coords array -------------------------------------------------------------------------- */ - -double DumpXTC::memory_usage() -{ - double bytes = maxbuf * sizeof(double); - bytes += 3*natoms * sizeof(float); - return bytes; + if (nevery_save == 0) nevery_save = output->every_dump[idump]; + else if (nevery_save != output->every_dump[idump]) + error->all("Cannot change dump_modify every for dump dcd"); } /* ---------------------------------------------------------------------- */ @@ -151,13 +136,11 @@ void DumpXTC::openfile() void DumpXTC::write_header(int n) { - // all procs realloc types & coords if necessary + // all procs realloc coords if total count grew if (n != natoms) { - memory->sfree(coords); - if (atom->tag_consecutive() == 0) - error->all("Atom IDs must be consecutive for dump xtc"); natoms = n; + memory->sfree(coords); coords = (float *) memory->smalloc(3*natoms*sizeof(float),"dump:coords"); } @@ -201,21 +184,32 @@ void DumpXTC::write_header(int n) int DumpXTC::count() { - return atom->nlocal; + if (igroup == 0) return atom->nlocal; + + int *mask = atom->mask; + int nlocal = atom->nlocal; + + int m = 0; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) m++; + return m; } /* ---------------------------------------------------------------------- */ -int DumpXTC::pack() +void DumpXTC::pack(int *ids) { + int m,n; + int *tag = atom->tag; double **x = atom->x; int *image = atom->image; + int *mask = atom->mask; int nlocal = atom->nlocal; // assume group all, so no need to perform mask check - int m = 0; + m = n = 0; if (unwrap_flag == 1) { double xprd = domain->xprd; double yprd = domain->yprd; @@ -224,55 +218,57 @@ int DumpXTC::pack() double xz = domain->xz; double yz = domain->yz; - for (int i = 0; i < nlocal; i++) { - int ix = (image[i] & 1023) - 512; - int iy = (image[i] >> 10 & 1023) - 512; - int iz = (image[i] >> 20) - 512; - - if (domain->triclinic) { - buf[m++] = tag[i]; - buf[m++] = sfactor * (x[i][0] + ix * xprd + iy * xy + iz * xz); - buf[m++] = sfactor * (x[i][1] + iy * yprd + iz * yz); - buf[m++] = sfactor * (x[i][2] + iz * zprd); - } else { - buf[m++] = tag[i]; - buf[m++] = sfactor * (x[i][0] + ix * xprd); - buf[m++] = sfactor * (x[i][1] + iy * yprd); - buf[m++] = sfactor * (x[i][2] + iz * zprd); + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + int ix = (image[i] & 1023) - 512; + int iy = (image[i] >> 10 & 1023) - 512; + int iz = (image[i] >> 20) - 512; + + if (domain->triclinic) { + buf[m++] = tag[i]; + buf[m++] = sfactor * (x[i][0] + ix * xprd + iy * xy + iz * xz); + buf[m++] = sfactor * (x[i][1] + iy * yprd + iz * yz); + buf[m++] = sfactor * (x[i][2] + iz * zprd); + } else { + buf[m++] = tag[i]; + buf[m++] = sfactor * (x[i][0] + ix * xprd); + buf[m++] = sfactor * (x[i][1] + iy * yprd); + buf[m++] = sfactor * (x[i][2] + iz * zprd); + } + ids[n++] = tag[i]; } - } - } else { - for (int i = 0; i < nlocal; i++) { - buf[m++] = tag[i]; - buf[m++] = sfactor*x[i][0]; - buf[m++] = sfactor*x[i][1]; - buf[m++] = sfactor*x[i][2]; - } - } - return m; + } else { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + buf[m++] = tag[i]; + buf[m++] = sfactor*x[i][0]; + buf[m++] = sfactor*x[i][1]; + buf[m++] = sfactor*x[i][2]; + ids[n++] = tag[i]; + } + } } /* ---------------------------------------------------------------------- */ void DumpXTC::write_data(int n, double *mybuf) { - float *xyz; - int j,tag; + int j; + + // copy buf atom coords into global array int m = 0; + int k = 3*ntotal; for (int i = 0; i < n; i++) { - tag = static_cast (mybuf[m]) - 1; - j = 3*tag; - coords[j++] = mybuf[m+1]; - coords[j++] = mybuf[m+2]; - coords[j] = mybuf[m+3]; - m += size_one; + coords[k++] = mybuf[m++]; + coords[k++] = mybuf[m++]; + coords[k++] = mybuf[m++]; + ntotal++; } // if last chunk of atoms in this snapshot, write global arrays to file - ntotal += n; if (ntotal == natoms) { write_frame(); ntotal = 0; @@ -281,6 +277,40 @@ void DumpXTC::write_data(int n, double *mybuf) /* ---------------------------------------------------------------------- */ +int DumpXTC::modify_param(int narg, char **arg) +{ + if (strcmp(arg[0],"unwrap") == 0) { + if (narg < 2) error->all("Illegal dump_modify command"); + if (strcmp(arg[1],"yes") == 0) unwrap_flag = 1; + else if (strcmp(arg[1],"no") == 0) unwrap_flag = 0; + else error->all("Illegal dump_modify command"); + return 2; + } else if (strcmp(arg[0],"precision") == 0) { + if (narg < 2) error->all("Illegal dump_modify command"); + precision = atof(arg[1]); + if ((fabs(precision-10.0) > EPS) && (fabs(precision-100.0) > EPS) && + (fabs(precision-1000.0) > EPS) && (fabs(precision-10000.0) > EPS) && + (fabs(precision-100000.0) > EPS) && + (fabs(precision-1000000.0) > EPS)) + error->all("Illegal dump_modify command"); + return 2; + } + return 0; +} + +/* ---------------------------------------------------------------------- + return # of bytes of allocated memory in buf and global coords array +------------------------------------------------------------------------- */ + +double DumpXTC::memory_usage() +{ + double bytes = Dump::memory_usage(); + bytes += 3*natoms * sizeof(float); + return bytes; +} + +/* ---------------------------------------------------------------------- */ + void DumpXTC::write_frame() { xdr3dfcoord(&xd,coords,&natoms,&precision); diff --git a/src/XTC/dump_xtc.h b/src/XTC/dump_xtc.h index 34392ea1d9..ebc928cc15 100644 --- a/src/XTC/dump_xtc.h +++ b/src/XTC/dump_xtc.h @@ -35,23 +35,24 @@ class DumpXTC : public Dump { public: DumpXTC(class LAMMPS *, int, char**); ~DumpXTC(); - void init(); - double memory_usage(); private: int natoms,ntotal; + int nevery_save; int unwrap_flag; // 1 if atom coords are unwrapped, 0 if no float precision; // user-adjustable precision setting float *coords; double sfactor; XDR xd; + void init_style(); int modify_param(int, char **); void openfile(); void write_header(int); int count(); - int pack(); + void pack(int *); void write_data(int, double *); + double memory_usage(); void write_frame(); }; diff --git a/src/XTC/xdr_compat.h b/src/XTC/xdr_compat.h index 37cad5e076..cb9be93991 100644 --- a/src/XTC/xdr_compat.h +++ b/src/XTC/xdr_compat.h @@ -60,6 +60,11 @@ extern "C" { typedef int bool_t; +#if defined(__MINGW32_VERSION) +typedef char * caddr_t; +typedef unsigned int u_int; +#endif + /* * Aninteger type that is 32 bits wide. Check if int, * long or short is 32 bits and die if none of them is :-) diff --git a/src/dump.cpp b/src/dump.cpp index 92ecb262d7..e740352921 100644 --- a/src/dump.cpp +++ b/src/dump.cpp @@ -17,6 +17,7 @@ #include "stdio.h" #include "dump.h" #include "atom.h" +#include "irregular.h" #include "update.h" #include "domain.h" #include "group.h" @@ -26,6 +27,17 @@ using namespace LAMMPS_NS; +// allocate space for static class variable + +Dump *Dump::dumpptr; + +#define BIG 1.0e20 +#define IBIG 2147483647 +#define EPSILON 1.0e-6 + +#define MIN(A,B) ((A) < (B)) ? (A) : (B) +#define MAX(A,B) ((A) > (B)) ? (A) : (B) + /* ---------------------------------------------------------------------- */ Dump::Dump(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) @@ -56,8 +68,10 @@ Dump::Dump(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) sort_flag = 0; append_flag = 0; - maxbuf = 0; - buf = NULL; + maxbuf = maxsort = maxproc = 0; + buf = bufsort = NULL; + idsort = index = proclist = NULL; + irregular = NULL; // parse filename for special syntax // if contains '%', write one file per proc and replace % with proc-ID @@ -108,6 +122,11 @@ Dump::~Dump() delete [] format_user; memory->sfree(buf); + memory->sfree(bufsort); + memory->sfree(idsort); + memory->sfree(index); + memory->sfree(proclist); + delete irregular; // XTC style sets fp to NULL since it closes file in its destructor @@ -124,6 +143,36 @@ Dump::~Dump() /* ---------------------------------------------------------------------- */ +void Dump::init() +{ + init_style(); + + if (!sort_flag) { + memory->sfree(buf); + memory->sfree(bufsort); + memory->sfree(idsort); + memory->sfree(index); + memory->sfree(proclist); + delete irregular; + + maxbuf = maxsort = maxproc = 0; + buf = bufsort = NULL; + idsort = index = proclist = NULL; + irregular = NULL; + } + + if (sort_flag && sortcol == 0 && atom->tag_enable == 0) + error->all("Cannot use dump sort on atom IDs with no atom IDs defined"); + + if (sort_flag && sortcol > size_one) + error->all("Dump sort column is invalid"); + + if (sort_flag && nprocs > 1 && irregular == NULL) + irregular = new Irregular(lmp); +} + +/* ---------------------------------------------------------------------- */ + void Dump::write() { // if file per timestep, open new file @@ -152,11 +201,12 @@ void Dump::write() } // nme = # of dump lines this proc will contribute to dump + + nme = count(); + // ntotal = total # of dump lines // nmax = max # of dump lines on any proc - int nme = count(); - int ntotal,nmax; if (multiproc) nmax = nme; else { @@ -169,25 +219,31 @@ void Dump::write() if (multiproc) write_header(nme); else write_header(ntotal); - // grow communication buffer if necessary + // pack my data into buf + // if sorting on IDs also request ID list from pack() + // sort buf as needed - if (nmax*size_one > maxbuf) { - maxbuf = nmax*size_one; + if (nmax > maxbuf) { + maxbuf = nmax; memory->sfree(buf); - buf = (double *) memory->smalloc(maxbuf*sizeof(double),"dump:buf"); + buf = (double *) + memory->smalloc(maxbuf*size_one*sizeof(double),"dump:buf"); + if (sort_flag && sortcol == 0) { + memory->sfree(ids); + ids = (int *) memory->smalloc(maxbuf*sizeof(int),"dump:ids"); + } } - // pack my data into buf - // me_size = # of quantities in buf - - int me_size = pack(); + if (sort_flag && sortcol == 0) pack(ids); + else pack(NULL); + if (sort_flag) sort(); // multiproc = 1 = each proc writes own data to own file // multiproc = 0 = all procs write to one file thru proc 0 // proc 0 pings each proc, receives it's data, writes to file // all other procs wait for ping, send their data to proc 0 - if (multiproc) write_data(me_size/size_one,buf); + if (multiproc) write_data(nme,buf); else { int tmp,nlines; MPI_Status status; @@ -201,7 +257,7 @@ void Dump::write() MPI_Wait(&request,&status); MPI_Get_count(&status,MPI_DOUBLE,&nlines); nlines /= size_one; - } else nlines = me_size/size_one; + } else nlines = nme; write_data(nlines,buf); } @@ -209,7 +265,7 @@ void Dump::write() } else { MPI_Recv(&tmp,0,MPI_INT,0,0,world,&status); - MPI_Rsend(buf,me_size,MPI_DOUBLE,0,0,world); + MPI_Rsend(buf,nme*size_one,MPI_DOUBLE,0,0,world); } } @@ -226,6 +282,229 @@ void Dump::write() } } +/* ---------------------------------------------------------------------- + generic opening of a dump file + ASCII or binary or gzipped + some derived classes override this function +------------------------------------------------------------------------- */ + +void Dump::openfile() +{ + // single file, already opened, so just return + + if (singlefile_opened) return; + if (multifile == 0) singlefile_opened = 1; + + // if one file per timestep, replace '*' with current timestep + + char *filecurrent; + if (multifile == 0) filecurrent = filename; + else { + filecurrent = new char[strlen(filename) + 16]; + char *ptr = strchr(filename,'*'); + *ptr = '\0'; + sprintf(filecurrent,"%s%d%s",filename,update->ntimestep,ptr+1); + *ptr = '*'; + } + + // open one file on proc 0 or file on every proc + + if (me == 0 || multiproc) { + if (compressed) { +#ifdef LAMMPS_GZIP + char gzip[128]; + sprintf(gzip,"gzip -6 > %s",filecurrent); + fp = popen(gzip,"w"); +#else + error->one("Cannot open gzipped file"); +#endif + } else if (binary) { + fp = fopen(filecurrent,"wb"); + } else if (append_flag) { + fp = fopen(filecurrent,"a"); + } else { + fp = fopen(filecurrent,"w"); + } + + if (fp == NULL) error->one("Cannot open dump file"); + } else fp = NULL; + + // delete string with timestep replaced + + if (multifile) delete [] filecurrent; +} + +/* ---------------------------------------------------------------------- + parallel sort of buf across all procs + changes nme, reorders datums in buf, grows buf if necessary +------------------------------------------------------------------------- */ + +void Dump::sort() +{ + int i,iproc; + double value; + + // if single proc, swap ptrs to buf,ids <-> bufsort,idsort + + if (nprocs == 1) { + if (nme > maxsort) { + maxsort = nme; + memory->sfree(bufsort); + bufsort = (double *) + memory->smalloc(maxsort*size_one*sizeof(double),"dump:bufsort"); + memory->sfree(index); + index = (int *) memory->smalloc(maxsort*sizeof(int),"dump:index"); + if (sortcol == 0) { + memory->sfree(idsort); + idsort = (int *) memory->smalloc(maxsort*sizeof(int),"dump:idsort"); + } + } + + double *dptr = buf; + buf = bufsort; + bufsort = dptr; + + if (sortcol == 0) { + int *iptr = ids; + ids = idsort; + idsort = iptr; + } + + // if multiple procs, exchange datums between procs via irregular + + } else { + + // grow proclist if necessary + + if (nme > maxproc) { + maxproc = nme; + memory->sfree(proclist); + proclist = (int *) memory->smalloc(maxproc*sizeof(int),"dump:proclist"); + } + + // proclist[i] = which proc Ith datum will be sent to + + if (sortcol == 0) { + int min = IBIG; + int max = 0; + for (i = 0; i < nme; i++) { + min = MIN(min,ids[i]); + max = MAX(max,ids[i]); + } + int minall,maxall; + MPI_Allreduce(&min,&minall,1,MPI_INT,MPI_MIN,world); + MPI_Allreduce(&max,&maxall,1,MPI_INT,MPI_MAX,world); + double range = maxall - minall + 0.1; + for (i = 0; i < nme; i++) { + iproc = static_cast ((ids[i]-minall)/range * nprocs); + proclist[i] = iproc; + } + + } else { + double min = BIG; + double max = -BIG; + for (i = 0; i < nme; i++) { + value = buf[i*size_one + sortcolm1]; + min = MIN(min,value); + max = MAX(max,value); + } + double minall,maxall; + MPI_Allreduce(&min,&minall,1,MPI_DOUBLE,MPI_MIN,world); + MPI_Allreduce(&max,&maxall,1,MPI_DOUBLE,MPI_MAX,world); + double range = maxall-minall + EPSILON*(maxall-minall); + if (range == 0.0) range = EPSILON; + for (i = 0; i < nme; i++) { + value = buf[i*size_one + sortcolm1]; + iproc = static_cast ((value-minall)/range * nprocs); + proclist[i] = iproc; + } + } + + // create comm plan, grow recv bufs if necessary, + // exchange datums, destroy plan + // if sorting on atom IDs, exchange IDs also + + nme = irregular->create_data(nme,proclist); + + if (nme > maxsort) { + maxsort = nme; + memory->sfree(bufsort); + bufsort = (double *) + memory->smalloc(maxsort*size_one*sizeof(double),"dump:bufsort"); + memory->sfree(index); + index = (int *) memory->smalloc(maxsort*sizeof(int),"dump:index"); + if (sortcol == 0) { + memory->sfree(idsort); + idsort = (int *) memory->smalloc(maxsort*sizeof(int),"dump:idsort"); + } + } + + irregular->exchange_data((char *) buf,size_one*sizeof(double), + (char *) bufsort); + if (sortcol == 0) + irregular->exchange_data((char *) ids,sizeof(int),(char *) idsort); + irregular->destroy_data(); + } + + // quicksort indices using IDs or buf column as comparator + + dumpptr = this; + for (i = 0; i < nme; i++) index[i] = i; + if (sortcol == 0) qsort(index,nme,sizeof(int),idcompare); + else qsort(index,nme,sizeof(int),bufcompare); + + // copy data from bufsort to buf using index + + if (nme > maxbuf) { + maxbuf = nme; + memory->sfree(buf); + buf = (double *) + memory->smalloc(maxbuf*size_one*sizeof(double),"dump:buf"); + } + + int nbytes = size_one*sizeof(double); + for (i = 0; i < nme; i++) + memcpy(&buf[i*size_one],&bufsort[index[i]*size_one],nbytes); +} + +/* ---------------------------------------------------------------------- + compare two atom IDs in sort data structure + called via qsort_r in sort() method + is a static method so access sort data structure via ptr +------------------------------------------------------------------------- */ + +int Dump::idcompare(const void *pi, const void *pj) +{ + int *idsort = dumpptr->idsort; + + int i = *((int *) pi); + int j = *((int *) pj); + + if (idsort[i] < idsort[j]) return -1; + if (idsort[i] > idsort[j]) return 1; + return 0; +} + +/* ---------------------------------------------------------------------- + compare two buffer quantities in sort data structure with size_one stride + called via qsort_r in sort() method + is a static method so access sort data structure via ptr +------------------------------------------------------------------------- */ + +int Dump::bufcompare(const void *pi, const void *pj) +{ + double *bufsort = dumpptr->bufsort; + int size_one = dumpptr->size_one; + int sortcolm1 = dumpptr->sortcolm1; + + int i = *((int *) pi)*size_one + sortcolm1; + int j = *((int *) pj)*size_one + sortcolm1; + + if (bufsort[i] < bufsort[j]) return -1; + if (bufsort[i] > bufsort[j]) return 1; + return 0; +} + /* ---------------------------------------------------------------------- process params common to all dumps here if unknown param, call modify_param specific to the dump @@ -285,9 +564,13 @@ void Dump::modify_params(int narg, char **arg) iarg += 2; } else if (strcmp(arg[iarg],"sort") == 0) { if (iarg+2 > narg) error->all("Illegal dump_modify command"); - if (strcmp(arg[iarg+1],"yes") == 0) sort_flag = 1; - else if (strcmp(arg[iarg+1],"no") == 0) sort_flag = 0; - else error->all("Illegal dump_modify command"); + if (strcmp(arg[iarg+1],"off") == 0) sort_flag = 0; + else { + sort_flag = 1; + sortcol = atoi(arg[iarg+1]); + if (sortcol < 0) error->all("Illegal dump_modify command"); + sortcolm1 = sortcol - 1; + } iarg += 2; } else { int n = modify_param(narg-iarg,&arg[iarg]); @@ -298,63 +581,19 @@ void Dump::modify_params(int narg, char **arg) } /* ---------------------------------------------------------------------- - return # of bytes of allocated memory in buf + return # of bytes of allocated memory ------------------------------------------------------------------------- */ double Dump::memory_usage() { - double bytes = maxbuf * sizeof(double); + double bytes = maxbuf*size_one * sizeof(double); + if (sort_flag) { + if (sortcol == 0) bytes += maxbuf * sizeof(int); // ids + bytes += maxsort*size_one * sizeof(double); // bufsort + bytes += maxsort * sizeof(int); // index + if (sortcol == 0) bytes += maxsort * sizeof(int); // idsort + bytes += maxproc * sizeof(int); // proclist + if (irregular) bytes += irregular->memory_usage(); + } return bytes; } - -/* ---------------------------------------------------------------------- - generic opening of a dump file - ASCII or binary or gzipped - some derived classes override this function -------------------------------------------------------------------------- */ - -void Dump::openfile() -{ - // single file, already opened, so just return - - if (singlefile_opened) return; - if (multifile == 0) singlefile_opened = 1; - - // if one file per timestep, replace '*' with current timestep - - char *filecurrent; - if (multifile == 0) filecurrent = filename; - else { - filecurrent = new char[strlen(filename) + 16]; - char *ptr = strchr(filename,'*'); - *ptr = '\0'; - sprintf(filecurrent,"%s%d%s",filename,update->ntimestep,ptr+1); - *ptr = '*'; - } - - // open one file on proc 0 or file on every proc - - if (me == 0 || multiproc) { - if (compressed) { -#ifdef LAMMPS_GZIP - char gzip[128]; - sprintf(gzip,"gzip -6 > %s",filecurrent); - fp = popen(gzip,"w"); -#else - error->one("Cannot open gzipped file"); -#endif - } else if (binary) { - fp = fopen(filecurrent,"wb"); - } else if (append_flag) { - fp = fopen(filecurrent,"a"); - } else { - fp = fopen(filecurrent,"w"); - } - - if (fp == NULL) error->one("Cannot open dump file"); - } else fp = NULL; - - // delete string with timestep replaced - - if (multifile) delete [] filecurrent; -} diff --git a/src/dump.h b/src/dump.h index b6029518bd..6e858ad9ca 100644 --- a/src/dump.h +++ b/src/dump.h @@ -24,6 +24,22 @@ class Dump : protected Pointers { char *id; // user-defined name of Dump char *style; // style of Dump int igroup,groupbit; // group that Dump is performed on + + int first_flag; // 0 if no initial dump, 1 if yes initial dump + int clearstep; // 1 if dump invokes computes, 0 if not + + // static variable across all Dump objects + + static Dump *dumpptr; // holds a ptr to Dump currently being used + + Dump(class LAMMPS *, int, char **); + virtual ~Dump(); + void init(); + void write(); + void modify_params(int, char **); + virtual double memory_usage(); + + protected: int me,nprocs; // proc info char *filename; // user-specified file @@ -33,41 +49,46 @@ class Dump : protected Pointers { int multiproc; // 0 = proc 0 writes for all, 1 = one file/proc int header_flag; // 0 = item, 2 = xyz - int first_flag; // 0 if no initial dump, 1 if yes initial dump int flush_flag; // 0 if no flush, 1 if flush every dump - int sort_flag; // 1 if write in sorted order, 0 if not + int sort_flag; // 1 if sorted output int append_flag; // 1 if open file in append mode, 0 if not int singlefile_opened; // 1 = one big file, already opened, else 0 - int clearstep; // 1 if dump invokes computes, 0 if not + int sortcol; // 0 to sort on ID, 1-N on columns + int sortcolm1; // sortcol - 1 char *format_default; // default format string char *format_user; // format string set by user char *format; // format string for the file write - double *buf; // memory for atom quantities - int maxbuf; // size of buf FILE *fp; // file to write dump to int size_one; // # of quantities for one atom + int nme; // # of atoms in this dump from me - Dump(class LAMMPS *, int, char **); - virtual ~Dump(); - virtual void init() {} - void write(); - void modify_params(int, char **); - virtual double memory_usage(); - - protected: double boxxlo,boxxhi; // local copies of domain values double boxylo,boxyhi; // lo/hi are bounding box for triclinic double boxzlo,boxzhi; double boxxy,boxxz,boxyz; + int maxbuf; // size of buf and ids + int maxsort; // size of bufsort, idsort, index + int maxproc; // size of proclist + double *buf; // memory for atom quantities + int *ids; // list of atom IDs, if sorting on IDs + double *bufsort; + int *idsort,*index,*proclist; + + class Irregular *irregular; + + virtual void init_style() = 0; virtual void openfile(); virtual int modify_param(int, char **) {return 0;} - virtual void write_header(int) = 0; virtual int count() = 0; - virtual int pack() = 0; + virtual void pack(int *) = 0; virtual void write_data(int, double *) = 0; + + void sort(); + static int idcompare(const void *, const void *); + static int bufcompare(const void *, const void *); }; } diff --git a/src/dump_atom.cpp b/src/dump_atom.cpp index 85a1701ccb..930fc6e1d3 100644 --- a/src/dump_atom.cpp +++ b/src/dump_atom.cpp @@ -34,7 +34,7 @@ DumpAtom::DumpAtom(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg) /* ---------------------------------------------------------------------- */ -void DumpAtom::init() +void DumpAtom::init_style() { if (image_flag == 0) size_one = 5; else size_one = 8; @@ -146,9 +146,9 @@ int DumpAtom::count() /* ---------------------------------------------------------------------- */ -int DumpAtom::pack() +void DumpAtom::pack(int *ids) { - return (this->*pack_choice)(); + (this->*pack_choice)(ids); } /* ---------------------------------------------------------------------- */ @@ -233,8 +233,10 @@ void DumpAtom::header_item_triclinic(int ndump) /* ---------------------------------------------------------------------- */ -int DumpAtom::pack_scale_image() +void DumpAtom::pack_scale_image(int *ids) { + int m,n; + int *tag = atom->tag; int *type = atom->type; int *image = atom->image; @@ -246,7 +248,7 @@ int DumpAtom::pack_scale_image() double invyprd = 1.0/domain->yprd; double invzprd = 1.0/domain->zprd; - int m = 0; + m = n = 00; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { buf[m++] = tag[i]; @@ -257,15 +259,16 @@ int DumpAtom::pack_scale_image() buf[m++] = (image[i] & 1023) - 512; buf[m++] = (image[i] >> 10 & 1023) - 512; buf[m++] = (image[i] >> 20) - 512; + if (ids) ids[n++] = tag[i]; } - - return m; } /* ---------------------------------------------------------------------- */ -int DumpAtom::pack_scale_noimage() +void DumpAtom::pack_scale_noimage(int *ids) { + int m,n; + int *tag = atom->tag; int *type = atom->type; int *mask = atom->mask; @@ -276,7 +279,7 @@ int DumpAtom::pack_scale_noimage() double invyprd = 1.0/domain->yprd; double invzprd = 1.0/domain->zprd; - int m = 0; + m = n = 0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { buf[m++] = tag[i]; @@ -284,15 +287,16 @@ int DumpAtom::pack_scale_noimage() buf[m++] = (x[i][0] - boxxlo) * invxprd; buf[m++] = (x[i][1] - boxylo) * invyprd; buf[m++] = (x[i][2] - boxzlo) * invzprd; + if (ids) ids[n++] = tag[i]; } - - return m; } /* ---------------------------------------------------------------------- */ -int DumpAtom::pack_scale_image_triclinic() +void DumpAtom::pack_scale_image_triclinic(int *ids) { + int m,n; + int *tag = atom->tag; int *type = atom->type; int *image = atom->image; @@ -302,7 +306,7 @@ int DumpAtom::pack_scale_image_triclinic() double lamda[3]; - int m = 0; + m = n = 0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { buf[m++] = tag[i]; @@ -314,15 +318,16 @@ int DumpAtom::pack_scale_image_triclinic() buf[m++] = (image[i] & 1023) - 512; buf[m++] = (image[i] >> 10 & 1023) - 512; buf[m++] = (image[i] >> 20) - 512; + if (ids) ids[n++] = tag[i]; } - - return m; } /* ---------------------------------------------------------------------- */ -int DumpAtom::pack_scale_noimage_triclinic() +void DumpAtom::pack_scale_noimage_triclinic(int *ids) { + int m,n; + int *tag = atom->tag; int *type = atom->type; int *mask = atom->mask; @@ -331,7 +336,7 @@ int DumpAtom::pack_scale_noimage_triclinic() double lamda[3]; - int m = 0; + m = n = 0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { buf[m++] = tag[i]; @@ -340,15 +345,16 @@ int DumpAtom::pack_scale_noimage_triclinic() buf[m++] = lamda[0]; buf[m++] = lamda[1]; buf[m++] = lamda[2]; + if (ids) ids[n++] = tag[i]; } - - return m; } /* ---------------------------------------------------------------------- */ -int DumpAtom::pack_noscale_image() +void DumpAtom::pack_noscale_image(int *ids) { + int m,n; + int *tag = atom->tag; int *type = atom->type; int *image = atom->image; @@ -356,7 +362,7 @@ int DumpAtom::pack_noscale_image() double **x = atom->x; int nlocal = atom->nlocal; - int m = 0; + m = n = 0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { buf[m++] = tag[i]; @@ -367,22 +373,23 @@ int DumpAtom::pack_noscale_image() buf[m++] = (image[i] & 1023) - 512; buf[m++] = (image[i] >> 10 & 1023) - 512; buf[m++] = (image[i] >> 20) - 512; + if (ids) ids[n++] = tag[i]; } - - return m; } /* ---------------------------------------------------------------------- */ -int DumpAtom::pack_noscale_noimage() +void DumpAtom::pack_noscale_noimage(int *ids) { + int m,n; + int *tag = atom->tag; int *type = atom->type; int *mask = atom->mask; double **x = atom->x; int nlocal = atom->nlocal; - int m = 0; + m = n = 0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { buf[m++] = tag[i]; @@ -390,9 +397,8 @@ int DumpAtom::pack_noscale_noimage() buf[m++] = x[i][0]; buf[m++] = x[i][1]; buf[m++] = x[i][2]; + if (ids) ids[n++] = tag[i]; } - - return m; } /* ---------------------------------------------------------------------- */ diff --git a/src/dump_atom.h b/src/dump_atom.h index 3120fca844..14f8e75248 100644 --- a/src/dump_atom.h +++ b/src/dump_atom.h @@ -27,7 +27,6 @@ namespace LAMMPS_NS { class DumpAtom : public Dump { public: DumpAtom(LAMMPS *, int, char**); - void init(); private: int scale_flag; // 1 if atom coords are scaled, 0 if no @@ -35,10 +34,11 @@ class DumpAtom : public Dump { char *columns; // column labels + void init_style(); int modify_param(int, char **); void write_header(int); int count(); - int pack(); + void pack(int *); void write_data(int, double *); typedef void (DumpAtom::*FnPtrHeader)(int); @@ -48,14 +48,14 @@ class DumpAtom : public Dump { void header_item(int); void header_item_triclinic(int); - typedef int (DumpAtom::*FnPtrPack)(); + typedef void (DumpAtom::*FnPtrPack)(int *); FnPtrPack pack_choice; // ptr to pack functions - int pack_scale_image(); - int pack_scale_noimage(); - int pack_noscale_image(); - int pack_noscale_noimage(); - int pack_scale_image_triclinic(); - int pack_scale_noimage_triclinic(); + void pack_scale_image(int *); + void pack_scale_noimage(int *); + void pack_noscale_image(int *); + void pack_noscale_noimage(int *); + void pack_scale_image_triclinic(int *); + void pack_scale_noimage_triclinic(int *); typedef void (DumpAtom::*FnPtrData)(int, double *); FnPtrData write_choice; // ptr to write data functions diff --git a/src/dump_cfg.cpp b/src/dump_cfg.cpp index 43548f364c..b11a5169b1 100755 --- a/src/dump_cfg.cpp +++ b/src/dump_cfg.cpp @@ -49,9 +49,6 @@ DumpCFG::DumpCFG(LAMMPS *lmp, int narg, char **arg) : // arrays for data rearrangement - recvcounts = new int[nprocs]; - displs = new int[nprocs]; - tags = NULL; rbuf = NULL; nchosen = nlines = 0; @@ -103,10 +100,6 @@ DumpCFG::~DumpCFG() delete [] typenames; } - delete [] recvcounts; - delete [] displs; - - if (tags) memory->sfree(tags); if (rbuf) memory->destroy_2d_double_array(rbuf); if (auxname) { @@ -117,7 +110,7 @@ DumpCFG::~DumpCFG() /* ---------------------------------------------------------------------- */ -void DumpCFG::init() +void DumpCFG::init_style() { if (multifile == 0) error->all("Dump in CFG format requires one snapshot per file"); @@ -222,10 +215,6 @@ void DumpCFG::write_header(int n) if (rbuf) memory->destroy_2d_double_array(rbuf); rbuf = memory->create_2d_double_array(nchosen,size_one,"dump:rbuf"); } - - // create a sorted list of atom IDs on writing proc(s) if necessary - - if (sort_flag) create_sorted_tags(); } /* ---------------------------------------------------------------------- @@ -239,147 +228,40 @@ void DumpCFG::write_data(int n, double *mybuf) int tag_i,index; double *mass = atom->mass; - // if sort flag is off, transfer data in buf to rbuf directly + // transfer data from buf to rbuf // if write by proc 0, transfer chunk by chunk - if (sort_flag == 0) { - for (i = 0, m = 0; i < n; i++) { - for (j = 0; j < size_one; j++) - rbuf[nlines][j] = mybuf[m++]; - nlines++; - } + for (i = 0, m = 0; i < n; i++) { + for (j = 0; j < size_one; j++) + rbuf[nlines][j] = mybuf[m++]; + nlines++; + } - // write data lines in rbuf to file after transfer is done - - if (nlines == nchosen) { - for (itype = 1; itype <= ntypes; itype++) { - for (i = 0; i < nchosen; i++) - if (rbuf[i][1] == itype) break; - if (i < nchosen) { - fprintf(fp,"%g\n",mass[itype]); - fprintf(fp,"%s\n",typenames[itype]); - for (; i < nchosen; i++) { - if (rbuf[i][1] == itype) { - for (j = 2; j < size_one; j++) { - if (vtype[j] == INT) - fprintf(fp,vformat[j],static_cast (rbuf[i][j])); - else fprintf(fp,vformat[j],rbuf[i][j]); - } - fprintf(fp,"\n"); + // write data lines in rbuf to file after transfer is done + + if (nlines == nchosen) { + for (itype = 1; itype <= ntypes; itype++) { + for (i = 0; i < nchosen; i++) + if (rbuf[i][1] == itype) break; + if (i < nchosen) { + fprintf(fp,"%g\n",mass[itype]); + fprintf(fp,"%s\n",typenames[itype]); + for (; i < nchosen; i++) { + if (rbuf[i][1] == itype) { + for (j = 2; j < size_one; j++) { + if (vtype[j] == INT) + fprintf(fp,vformat[j],static_cast (rbuf[i][j])); + else fprintf(fp,vformat[j],rbuf[i][j]); } + fprintf(fp,"\n"); } } } - nlines = 0; - } - - // if sort flag is on, transfer data lines in buf to rbuf & - // rearrange them in a sorted order - // if write by proc 0, transfer chunk by chunk - - } else { - - // sort data lines: - // find the index of the atom ID of each data line within the - // sorted atom IDs array - // transfer the data line to rbuf with its location according - // to this index - - for (i = 0, m = 0; i < n; i++) { - tag_i = static_cast(mybuf[m]); - for (j = 0; j < nchosen; j++) { - if (tag_i == tags[j]) { - index = j; - break; - } - } - for (j = 0; j < size_one; j++) - rbuf[index][j] = mybuf[m++]; - } - - // write data lines in rbuf to file after transfer is done - - nlines += n; - if (nlines == nchosen) { - for (itype = 1; itype <= ntypes; itype++) { - for (i = 0; i < nchosen; i++) - if (rbuf[i][1] == itype) break; - if (i < nchosen) { - fprintf(fp,"%g\n",mass[itype]); - fprintf(fp,"%s\n",typenames[itype]); - for (; i < nchosen; i++) { - if (rbuf[i][1] == itype) { - for (j = 2; j < size_one; j++) { - if (vtype[j] == INT) - fprintf(fp,vformat[j],static_cast(rbuf[i][j])); - else fprintf(fp,vformat[j],rbuf[i][j]); - } - fprintf(fp,"\n"); - } - } - } - } - nlines = 0; } + nlines = 0; } } -/* ---------------------------------------------------------------------- - create a sorted list of atom IDs on writing proc(s) -------------------------------------------------------------------------- */ - -void DumpCFG::create_sorted_tags() -{ - int index; - int *mytags; - int *tag = atom->tag; - int nlocal = atom->nlocal; - - index = 0; - - // if write by multiproc, each proc has its own atom IDs collected - - if (multiproc) { - if (tags) memory->sfree(tags); - tags = (int *) memory->smalloc(nchosen*sizeof(int),"dump:tags"); - for (int i = 0; i < nlocal; i++) - if (choose[i]) tags[index++] = tag[i]; - - // if write by proc 0, gather atom IDs of all the chosen atoms on proc 0 - - } else { - mytags = (int *) memory->smalloc(nmine*sizeof(int),"dump:mytags"); - for (int i = 0; i < nlocal; i++) - if (choose[i]) mytags[index++] = tag[i]; - MPI_Gather(&nmine,1,MPI_INT,recvcounts,1,MPI_INT,0,world); - if (me == 0) { - if (tags) memory->sfree(tags); - tags = (int *) memory->smalloc(nchosen*sizeof(int),"dump:tags"); - displs[0] = 0; - for (int iproc = 1; iproc < nprocs; iproc++) - displs[iproc] = displs[iproc-1] + recvcounts[iproc-1]; - } - MPI_Gatherv(mytags,nmine,MPI_INT,tags,recvcounts,displs,MPI_INT,0,world); - memory->sfree(mytags); - } - - if (multiproc || me == 0) qsort(tags,nchosen,sizeof(int),tag_compare); -} - -/* ---------------------------------------------------------------------- - compare function for qsort -------------------------------------------------------------------------- */ - -int DumpCFG::tag_compare(const void *itag, const void *jtag) -{ - int tag_i = *((int *) itag); - int tag_j = *((int *) jtag); - - if (tag_i < tag_j) return -1; - else if (tag_i > tag_j) return 1; - return 0; -} - /* ---------------------------------------------------------------------- */ int DumpCFG::modify_param2(int narg, char **arg) diff --git a/src/dump_cfg.h b/src/dump_cfg.h index 3abdfaa0ad..ca3cd4032c 100755 --- a/src/dump_cfg.h +++ b/src/dump_cfg.h @@ -28,10 +28,6 @@ class DumpCFG : public DumpCustom { public: DumpCFG(class LAMMPS *, int, char **); ~DumpCFG(); - void init(); - void write_header(int); - void write_data(int, double *); - int modify_param2(int, char **); private: int ntypes; // # of atom types @@ -39,13 +35,13 @@ class DumpCFG : public DumpCustom { char **auxname; // name strings of auxiliary properties int nchosen; // # of lines to be written on a writing proc int nlines; // # of lines transferred from buf to rbuf - int *recvcounts; // # of data (tag info) received from all procs - int *displs; // displacement of data (tag info) within buffer - int *tags; // atom IDs of all the chosen atoms collected double **rbuf; // buf of data lines for data lines rearrangement - void create_sorted_tags(); - static int tag_compare(const void *, const void *); + void init_style(); + void write_header(int); + void write_data(int, double *); + + int modify_param2(int, char **); }; } diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp index b385da6577..c5058d89fd 100644 --- a/src/dump_custom.cpp +++ b/src/dump_custom.cpp @@ -159,7 +159,7 @@ DumpCustom::~DumpCustom() /* ---------------------------------------------------------------------- */ -void DumpCustom::init() +void DumpCustom::init_style() { delete [] format; char *str; @@ -779,10 +779,16 @@ int DumpCustom::count() /* ---------------------------------------------------------------------- */ -int DumpCustom::pack() +void DumpCustom::pack(int *ids) { for (int n = 0; n < size_one; n++) (this->*pack_choice[n])(n); - return nmine*size_one; + if (ids) { + int *tag = atom->tag; + int nlocal = atom->nlocal; + int n = 0; + for (int i = 0; i < nlocal; i++) + if (choose[i]) ids[n++] = tag[i]; + } } /* ---------------------------------------------------------------------- */ @@ -1437,7 +1443,7 @@ int DumpCustom::modify_param(int narg, char **arg) nthresh++; return 4; - // pass along params to child class + // pass along unknown params to child class } else { int n = modify_param2(narg,arg); @@ -1453,7 +1459,7 @@ int DumpCustom::modify_param(int narg, char **arg) double DumpCustom::memory_usage() { - double bytes = maxbuf * sizeof(double); + double bytes = Dump::memory_usage(); bytes += maxlocal * sizeof(int); bytes += maxlocal * sizeof(double); bytes += maxlocal * nvariable * sizeof(double); diff --git a/src/dump_custom.h b/src/dump_custom.h index 158b52699b..18be848f00 100644 --- a/src/dump_custom.h +++ b/src/dump_custom.h @@ -28,8 +28,6 @@ class DumpCustom : public Dump { public: DumpCustom(class LAMMPS *, int, char **); virtual ~DumpCustom(); - virtual void init(); - double memory_usage(); protected: int nevery; // dump frequency to check Fix against @@ -71,10 +69,12 @@ class DumpCustom : public Dump { // private methods + virtual void init_style(); virtual void write_header(int); int count(); - int pack(); + void pack(int *); virtual void write_data(int, double *); + double memory_usage(); void parse_fields(int, char **); int add_compute(char *); diff --git a/src/dump_dcd.cpp b/src/dump_dcd.cpp index 77bf91c4dc..c0f4bc06c8 100644 --- a/src/dump_dcd.cpp +++ b/src/dump_dcd.cpp @@ -13,6 +13,7 @@ /* ---------------------------------------------------------------------- Contributing author: Naveen Michaud-Agrawal (Johns Hopkins U) + Axel Kohlmeyer (Temple U), support for groups ------------------------------------------------------------------------- */ #include "math.h" @@ -54,20 +55,22 @@ static inline void fwrite_int32(FILE* fd, uint32_t i) DumpDCD::DumpDCD(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg) { if (narg != 5) error->all("Illegal dump dcd command"); - if (igroup != group->find("all")) error->all("Dump dcd must use group all"); if (binary || compressed || multifile || multiproc) error->all("Invalid dump dcd filename"); - size_one = 4; + size_one = 3; + sort_flag = 1; + sortcol = 0; + unwrap_flag = 0; format_default = NULL; // allocate global array for atom coords - natoms = static_cast (atom->natoms); + if (igroup == 0) natoms = static_cast (atom->natoms); + else natoms = static_cast (group->count(igroup)); if (natoms <= 0) error->all("Invalid natoms for dump dcd"); - if (atom->tag_consecutive() == 0) - error->all("Atom IDs must be consecutive for dump dcd"); + coords = (float *) memory->smalloc(3*natoms*sizeof(float),"dump:coords"); xf = &coords[0*natoms]; yf = &coords[1*natoms]; @@ -88,8 +91,11 @@ DumpDCD::~DumpDCD() /* ---------------------------------------------------------------------- */ -void DumpDCD::init() +void DumpDCD::init_style() { + if (sort_flag == 0 || sortcol != 0) + error->all("Dump dcd requires sorting by atom ID"); + // check that dump frequency has not changed and is not a variable int idump; @@ -105,31 +111,6 @@ void DumpDCD::init() /* ---------------------------------------------------------------------- */ -int DumpDCD::modify_param(int narg, char **arg) -{ - if (strcmp(arg[0],"unwrap") == 0) { - if (narg < 2) error->all("Illegal dump_modify command"); - if (strcmp(arg[1],"yes") == 0) unwrap_flag = 1; - else if (strcmp(arg[1],"no") == 0) unwrap_flag = 0; - else error->all("Illegal dump_modify command"); - return 2; - } - return 0; -} - -/* ---------------------------------------------------------------------- - return # of bytes of allocated memory in buf and global coords array -------------------------------------------------------------------------- */ - -double DumpDCD::memory_usage() -{ - double bytes = maxbuf * sizeof(double); - bytes += 3*natoms * sizeof(float); - return bytes; -} - -/* ---------------------------------------------------------------------- */ - void DumpDCD::openfile() { if (me == 0) { @@ -193,21 +174,30 @@ void DumpDCD::write_header(int n) int DumpDCD::count() { - return atom->nlocal; + if (igroup == 0) return atom->nlocal; + + int *mask = atom->mask; + int nlocal = atom->nlocal; + + int m = 0; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) m++; + return m; } /* ---------------------------------------------------------------------- */ -int DumpDCD::pack() +void DumpDCD::pack(int *ids) { + int m,n; + int *tag = atom->tag; double **x = atom->x; int *image = atom->image; + int *mask = atom->mask; int nlocal = atom->nlocal; - // assume group all, so no need to perform mask check - - int m = 0; + m = n = 0; if (unwrap_flag) { double xprd = domain->xprd; double yprd = domain->yprd; @@ -217,53 +207,51 @@ int DumpDCD::pack() double yz = domain->yz; for (int i = 0; i < nlocal; i++) { - int ix = (image[i] & 1023) - 512; - int iy = (image[i] >> 10 & 1023) - 512; - int iz = (image[i] >> 20) - 512; + if (mask[i] & groupbit) { + int ix = (image[i] & 1023) - 512; + int iy = (image[i] >> 10 & 1023) - 512; + int iz = (image[i] >> 20) - 512; - if (domain->triclinic) { - buf[m++] = tag[i]; - buf[m++] = x[i][0] + ix * xprd + iy * xy + iz * xz; - buf[m++] = x[i][1] + iy * yprd + iz * yz; - buf[m++] = x[i][2] + iz * zprd; - } else { - buf[m++] = tag[i]; - buf[m++] = x[i][0] + ix * xprd; - buf[m++] = x[i][1] + iy * yprd; - buf[m++] = x[i][2] + iz * zprd; + if (domain->triclinic) { + buf[m++] = x[i][0] + ix * xprd + iy * xy + iz * xz; + buf[m++] = x[i][1] + iy * yprd + iz * yz; + buf[m++] = x[i][2] + iz * zprd; + } else { + buf[m++] = x[i][0] + ix * xprd; + buf[m++] = x[i][1] + iy * yprd; + buf[m++] = x[i][2] + iz * zprd; + } + ids[n++] = tag[i]; } } - } else { - for (int i = 0; i < nlocal; i++) { - buf[m++] = tag[i]; - buf[m++] = x[i][0]; - buf[m++] = x[i][1]; - buf[m++] = x[i][2]; - } - } - return m; + } else { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + buf[m++] = x[i][0]; + buf[m++] = x[i][1]; + buf[m++] = x[i][2]; + ids[n++] = tag[i]; + } + } } /* ---------------------------------------------------------------------- */ void DumpDCD::write_data(int n, double *mybuf) { - // spread buf atom coords into global arrays + // copy buf atom coords into 3 global arrays - int tag; int m = 0; for (int i = 0; i < n; i++) { - tag = static_cast (mybuf[m]) - 1; - xf[tag] = mybuf[m+1]; - yf[tag] = mybuf[m+2]; - zf[tag] = mybuf[m+3]; - m += size_one; + xf[ntotal] = mybuf[m++]; + yf[ntotal] = mybuf[m++]; + zf[ntotal] = mybuf[m++]; + ntotal++; } // if last chunk of atoms in this snapshot, write global arrays to file - ntotal += n; if (ntotal == natoms) { write_frame(); ntotal = 0; @@ -272,6 +260,31 @@ void DumpDCD::write_data(int n, double *mybuf) /* ---------------------------------------------------------------------- */ +int DumpDCD::modify_param(int narg, char **arg) +{ + if (strcmp(arg[0],"unwrap") == 0) { + if (narg < 2) error->all("Illegal dump_modify command"); + if (strcmp(arg[1],"yes") == 0) unwrap_flag = 1; + else if (strcmp(arg[1],"no") == 0) unwrap_flag = 0; + else error->all("Illegal dump_modify command"); + return 2; + } + return 0; +} + +/* ---------------------------------------------------------------------- + return # of bytes of allocated memory in buf and global coords array +------------------------------------------------------------------------- */ + +double DumpDCD::memory_usage() +{ + double bytes = Dump::memory_usage(); + bytes += 3*natoms * sizeof(float); + return bytes; +} + +/* ---------------------------------------------------------------------- */ + void DumpDCD::write_frame() { // write coords diff --git a/src/dump_dcd.h b/src/dump_dcd.h index dc73483623..c92edc4436 100644 --- a/src/dump_dcd.h +++ b/src/dump_dcd.h @@ -30,21 +30,21 @@ class DumpDCD : public Dump { public: DumpDCD(LAMMPS *, int, char**); ~DumpDCD(); - void init(); - double memory_usage(); private: int natoms,ntotal,headerflag,nevery_save,nframes; float *coords,*xf,*yf,*zf; int unwrap_flag; // 1 if atom coords are unwrapped, 0 if no + void init_style(); void openfile(); void write_header(int); int count(); - int pack(); + void pack(int *); void write_data(int, double *); - int modify_param(int, char **); + double memory_usage(); + void write_frame(); void write_dcd_header(const char *); }; diff --git a/src/dump_local.cpp b/src/dump_local.cpp index bda5c599af..c93beac713 100644 --- a/src/dump_local.cpp +++ b/src/dump_local.cpp @@ -120,8 +120,11 @@ DumpLocal::~DumpLocal() /* ---------------------------------------------------------------------- */ -void DumpLocal::init() +void DumpLocal::init_style() { + if (sort_flag && sortcol == 0) + error->all("Dump local cannot sort by atom ID"); + delete [] format; char *str; if (format_user) str = format_user; @@ -237,10 +240,9 @@ int DumpLocal::count() /* ---------------------------------------------------------------------- */ -int DumpLocal::pack() +void DumpLocal::pack(int *dummy) { for (int n = 0; n < size_one; n++) (this->*pack_choice[n])(n); - return nmine*size_one; } /* ---------------------------------------------------------------------- */ diff --git a/src/dump_local.h b/src/dump_local.h index 43f2201766..53173d7910 100644 --- a/src/dump_local.h +++ b/src/dump_local.h @@ -28,7 +28,6 @@ class DumpLocal : public Dump { public: DumpLocal(LAMMPS *, int, char **); ~DumpLocal(); - void init(); private: int nevery; // dump frequency to check Fix against @@ -54,10 +53,11 @@ class DumpLocal : public Dump { char **id_fix; // their IDs class Fix **fix; // list of ptrs to the Fix objects + void init_style(); int modify_param(int, char **); void write_header(int); int count(); - int pack(); + void pack(int *); void write_data(int, double *); void parse_fields(int, char **); diff --git a/src/dump_xyz.cpp b/src/dump_xyz.cpp index 9219451fb3..93d6294aa8 100644 --- a/src/dump_xyz.cpp +++ b/src/dump_xyz.cpp @@ -28,39 +28,18 @@ DumpXYZ::DumpXYZ(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg) if (binary || multiproc) error->all("Invalid dump xyz filename"); size_one = 5; + sort_flag = 1; + sortcol = 0; char *str = (char *) "%d %g %g %g"; int n = strlen(str) + 1; format_default = new char[n]; strcpy(format_default,str); - - // allocate global array for atom coords if group is all - - if (igroup == 0) { - natoms = static_cast (atom->natoms); - if (natoms <= 0) error->all("Invalid natoms for dump xyz"); - if (atom->tag_consecutive() == 0) - error->all("Atom IDs must be consecutive for dump xyz"); - types = (int *) memory->smalloc(natoms*sizeof(int),"dump:types"); - coords = (float *) memory->smalloc(3*natoms*sizeof(float),"dump:coords"); - } - - ntotal = 0; } /* ---------------------------------------------------------------------- */ -DumpXYZ::~DumpXYZ() -{ - if (igroup == 0) { - memory->sfree(types); - memory->sfree(coords); - } -} - -/* ---------------------------------------------------------------------- */ - -void DumpXYZ::init() +void DumpXYZ::init_style() { delete [] format; char *str; @@ -77,41 +56,10 @@ void DumpXYZ::init() if (multifile == 0) openfile(); } -/* ---------------------------------------------------------------------- - return # of bytes of allocated memory in buf and global coords array -------------------------------------------------------------------------- */ - -double DumpXYZ::memory_usage() -{ - double bytes = maxbuf * sizeof(double); - if (igroup == 0) { - bytes += natoms * sizeof(int); - bytes += 3*natoms * sizeof(float); - } - return bytes; -} - /* ---------------------------------------------------------------------- */ void DumpXYZ::write_header(int n) { - // check bounds of atom IDs if group is all - // all procs realloc types & coords if necessary - - if (igroup == 0) { - if (atom->tag_consecutive() == 0) - error->all("Atom IDs must be consecutive for dump xyz all"); - if (n != natoms) { - memory->sfree(types); - memory->sfree(coords); - natoms = n; - types = (int *) memory->smalloc(natoms*sizeof(int),"dump:types"); - coords = (float *) memory->smalloc(3*natoms*sizeof(float),"dump:coords"); - } - } - - // only proc 0 writes header - if (me == 0) { fprintf(fp,"%d\n",n); fprintf(fp,"Atoms\n"); @@ -135,15 +83,17 @@ int DumpXYZ::count() /* ---------------------------------------------------------------------- */ -int DumpXYZ::pack() +void DumpXYZ::pack(int *ids) { + int m,n; + int *tag = atom->tag; int *type = atom->type; int *mask = atom->mask; double **x = atom->x; int nlocal = atom->nlocal; - int m = 0; + m = n = 0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { buf[m++] = tag[i]; @@ -151,57 +101,18 @@ int DumpXYZ::pack() buf[m++] = x[i][0]; buf[m++] = x[i][1]; buf[m++] = x[i][2]; + if (ids) ids[n++] = tag[i]; } - - return m; } /* ---------------------------------------------------------------------- */ void DumpXYZ::write_data(int n, double *mybuf) -{ - // for group = all, spread buf atom coords into global arrays - // if last chunk of atoms in this snapshot, write global arrays to file - - if (igroup == 0) { - int j,tag; - - int m = 0; - for (int i = 0; i < n; i++) { - tag = static_cast (mybuf[m]) - 1; - types[tag] = static_cast (mybuf[m+1]); - j = 3*tag; - coords[j++] = mybuf[m+2]; - coords[j++] = mybuf[m+3]; - coords[j] = mybuf[m+4]; - m += size_one; - } - - ntotal += n; - if (ntotal == natoms) { - write_frame(); - ntotal = 0; - } - - // for group != all, write unsorted type,x,y,z to file - - } else { - int m = 0; - for (int i = 0; i < n; i++) { - fprintf(fp,format, - static_cast (mybuf[m+1]),mybuf[m+2],mybuf[m+3],mybuf[m+4]); - m += size_one; - } - } -} - -/* ---------------------------------------------------------------------- */ - -void DumpXYZ::write_frame() { int m = 0; - for (int i = 0; i < natoms; i++) { - fprintf(fp,format,types[i],coords[m],coords[m+1],coords[m+2]); - m += 3; + for (int i = 0; i < n; i++) { + fprintf(fp,format, + static_cast (mybuf[m+1]),mybuf[m+2],mybuf[m+3],mybuf[m+4]); + m += size_one; } } diff --git a/src/dump_xyz.h b/src/dump_xyz.h index dbf4ff96ba..b155805113 100644 --- a/src/dump_xyz.h +++ b/src/dump_xyz.h @@ -27,20 +27,14 @@ namespace LAMMPS_NS { class DumpXYZ : public Dump { public: DumpXYZ(class LAMMPS *, int, char**); - ~DumpXYZ(); - void init(); - double memory_usage(); + ~DumpXYZ() {} private: - int natoms,ntotal; - int *types; - float *coords; - + void init_style(); void write_header(int); int count(); - int pack(); + void pack(int *); void write_data(int, double *); - void write_frame(); }; } diff --git a/src/irregular.cpp b/src/irregular.cpp index ded20b4efc..37af76759c 100644 --- a/src/irregular.cpp +++ b/src/irregular.cpp @@ -42,7 +42,11 @@ Irregular::Irregular(LAMMPS *lmp) : Pointers(lmp) procgrid = comm->procgrid; grid2proc = comm->grid2proc; - // initialize comm buffers + aplan = NULL; + dplan = NULL; + + // initialize buffers for atom comm, not used for datum comm + // these can persist for multiple irregular operations maxsend = BUFMIN; buf_send = (double *) @@ -56,6 +60,9 @@ Irregular::Irregular(LAMMPS *lmp) : Pointers(lmp) Irregular::~Irregular() { + if (aplan) destroy_atom(); + if (dplan) destroy_data(); + memory->sfree(buf_send); memory->sfree(buf_recv); } @@ -64,7 +71,7 @@ Irregular::~Irregular() communicate atoms to new owning procs via irregular communication can be used in place of comm->exchange() unlike exchange(), allows atoms to have moved arbitrarily long distances - first setup irregular comm pattern, then invoke it + sets up irregular plan, invokes it, destroys it atoms must be remapped to be inside simulation box before this is called for triclinic: atoms must be in lamda coords (0-1) before this is called ------------------------------------------------------------------------- */ @@ -127,11 +134,10 @@ void Irregular::migrate_atoms() // create irregular communication plan, perform comm, destroy plan // returned nrecv = size of buffer needed for incoming atoms - int nrecv; - PlanAtom *plan = create_atom(nsendatom,sizes,proclist,&nrecv); + int nrecv = create_atom(nsendatom,sizes,proclist); if (nrecv > maxrecv) grow_recv(nrecv); - exchange_atom(plan,buf_send,sizes,buf_recv); - destroy_atom(plan); + exchange_atom(buf_send,sizes,buf_recv); + destroy_atom(); delete [] sizes; delete [] proclist; @@ -147,22 +153,22 @@ void Irregular::migrate_atoms() } /* ---------------------------------------------------------------------- - create an irregular communication plan for atoms + create a communication plan for atoms n = # of atoms to send sizes = # of doubles for each atom - proclist = proc to send each atom to (none to me) - return nrecvsize = total # of doubles I will recv + proclist = proc to send each atom to (not including self) + return total # of doubles I will recv (not including self) ------------------------------------------------------------------------- */ -Irregular::PlanAtom *Irregular::create_atom(int n, int *sizes, - int *proclist, int *nrecvsize) +int Irregular::create_atom(int n, int *sizes, int *proclist) { int i; // allocate plan and work vectors - PlanAtom *plan = (struct PlanAtom *) - memory->smalloc(sizeof(PlanAtom),"irregular:plan"); + if (aplan) destroy_atom(); + aplan = (struct PlanAtom *) + memory->smalloc(sizeof(PlanAtom),"irregular:aplan"); int *list = new int[nprocs]; int *count = new int[nprocs]; @@ -221,7 +227,7 @@ Irregular::PlanAtom *Irregular::create_atom(int n, int *sizes, } } - // num_send = # of datums I send to each proc + // num_send = # of atoms I send to each proc for (i = 0; i < nsend; i++) num_send[i] = 0; for (i = 0; i < n; i++) { @@ -229,11 +235,11 @@ Irregular::PlanAtom *Irregular::create_atom(int n, int *sizes, num_send[isend]++; } - // count = offsets into n-length index_send for each proc I send to - // index_send = list of which datums to send to each proc - // 1st N1 values are datum indices for 1st proc, - // next N2 values are datum indices for 2nd proc, etc - // offset_send = where each datum starts in send buffer + // count = offsets into index_send for each proc I send to + // index_send = list of which atoms to send to each proc + // 1st N1 values are atom indices for 1st proc, + // next N2 values are atom indices for 2nd proc, etc + // offset_send = where each atom starts in send buffer count[0] = 0; for (i = 1; i < nsend; i++) count[i] = count[i-1] + num_send[i-1]; @@ -246,7 +252,7 @@ Irregular::PlanAtom *Irregular::create_atom(int n, int *sizes, } // tell receivers how much data I send - // sendmax = largest # of datums I send in a single message + // sendmax = largest # of doubles I send in a single message int sendmax = 0; for (i = 0; i < nsend; i++) { @@ -259,15 +265,15 @@ Irregular::PlanAtom *Irregular::create_atom(int n, int *sizes, // length_recv = total size of message each proc sends me // nrecvsize = total size of data I recv - *nrecvsize = 0; + int nrecvsize = 0; for (i = 0; i < nrecv; i++) { MPI_Recv(&length_recv[i],1,MPI_INT,MPI_ANY_SOURCE,0,world,status); proc_recv[i] = status->MPI_SOURCE; - *nrecvsize += length_recv[i]; + nrecvsize += length_recv[i]; } // barrier to insure all MPI_ANY_SOURCE messages are received - // else another proc could proceed to irregular_perform() and send to me + // else another proc could proceed to exchange_atom() and send to me MPI_Barrier(world); @@ -276,41 +282,115 @@ Irregular::PlanAtom *Irregular::create_atom(int n, int *sizes, delete [] count; delete [] list; - // initialize plan and return it + // initialize plan - plan->nsend = nsend; - plan->nrecv = nrecv; - plan->sendmax = sendmax; + aplan->nsend = nsend; + aplan->nrecv = nrecv; + aplan->sendmax = sendmax; - plan->proc_send = proc_send; - plan->length_send = length_send; - plan->num_send = num_send; - plan->index_send = index_send; - plan->offset_send = offset_send; - plan->proc_recv = proc_recv; - plan->length_recv = length_recv; - plan->request = request; - plan->status = status; + aplan->proc_send = proc_send; + aplan->length_send = length_send; + aplan->num_send = num_send; + aplan->index_send = index_send; + aplan->offset_send = offset_send; + aplan->proc_recv = proc_recv; + aplan->length_recv = length_recv; - return plan; + aplan->request = request; + aplan->status = status; + + return nrecvsize; } /* ---------------------------------------------------------------------- - create an irregular communication plan for datums - n = # of datums to send - proclist = proc to send each datum to (none to me) - return nrecvsize = total # of datums I will recv + communicate atoms via PlanAtom + sendbuf = list of atoms to send + sizes = # of doubles for each atom + recvbuf = received atoms ------------------------------------------------------------------------- */ -Irregular::PlanData *Irregular::create_data(int n, int *proclist, - int *nrecvsize) +void Irregular::exchange_atom(double *sendbuf, int *sizes, double *recvbuf) { - int i; + int i,m,n,offset,num_send; + + // post all receives + + offset = 0; + for (int irecv = 0; irecv < aplan->nrecv; irecv++) { + MPI_Irecv(&recvbuf[offset],aplan->length_recv[irecv],MPI_DOUBLE, + aplan->proc_recv[irecv],0,world,&aplan->request[irecv]); + offset += aplan->length_recv[irecv]; + } + + // allocate buf for largest send + + double *buf = (double *) memory->smalloc(aplan->sendmax*sizeof(double), + "irregular:buf"); + + // send each message + // pack buf with list of atoms + // m = index of atom in sendbuf + + int *index_send = aplan->index_send; + int nsend = aplan->nsend; + n = 0; + + for (int isend = 0; isend < nsend; isend++) { + offset = 0; + num_send = aplan->num_send[isend]; + for (i = 0; i < num_send; i++) { + m = index_send[n++]; + memcpy(&buf[offset],&sendbuf[aplan->offset_send[m]], + sizes[m]*sizeof(double)); + offset += sizes[m]; + } + MPI_Send(buf,aplan->length_send[isend],MPI_DOUBLE, + aplan->proc_send[isend],0,world); + } + + // free temporary send buffer + + memory->sfree(buf); + + // wait on all incoming messages + + if (aplan->nrecv) MPI_Waitall(aplan->nrecv,aplan->request,aplan->status); +} + +/* ---------------------------------------------------------------------- + destroy communication plan for atoms +------------------------------------------------------------------------- */ + +void Irregular::destroy_atom() +{ + delete [] aplan->proc_send; + delete [] aplan->length_send; + delete [] aplan->num_send; + delete [] aplan->index_send; + delete [] aplan->offset_send; + delete [] aplan->proc_recv; + delete [] aplan->length_recv; + delete [] aplan->request; + delete [] aplan->status; + memory->sfree(aplan); + aplan = NULL; +} + +/* ---------------------------------------------------------------------- + create a communication plan for datums + n = # of datums to send + proclist = proc to send each datum to (including self) + return total # of datums I will recv (including self) +------------------------------------------------------------------------- */ + +int Irregular::create_data(int n, int *proclist) +{ + int i,m; // allocate plan and work vectors - PlanData *plan = (struct PlanData *) - memory->smalloc(sizeof(PlanData),"irregular:plan"); + dplan = (struct PlanData *) + memory->smalloc(sizeof(PlanData),"irregular:dplan"); int *list = new int[nprocs]; int *count = new int[nprocs]; @@ -324,6 +404,7 @@ Irregular::PlanData *Irregular::create_data(int n, int *proclist, int nrecv; MPI_Reduce_scatter(list,&nrecv,count,MPI_INT,MPI_SUM,world); + if (list[me]) nrecv--; // allocate receive arrays @@ -340,42 +421,56 @@ Irregular::PlanData *Irregular::create_data(int n, int *proclist, int nsend = 0; for (i = 0; i < nprocs; i++) if (list[i]) nsend++; + if (list[me]) nsend--; - // allocate send arrays + // allocate send and self arrays int *proc_send = new int[nsend]; int *num_send = new int[nsend]; - int *index_send = new int[n]; + int *index_send = new int[n-list[me]]; + int *index_self = new int[list[me]]; // proc_send = procs I send to // num_send = # of datums I send to each proc + // num_self = # of datums I copy to self // to balance pattern of send messages: // each proc begins with iproc > me, continues until iproc = me + // reset list to store which send message each proc corresponds to + + int num_self; int iproc = me; int isend = 0; for (i = 0; i < nprocs; i++) { iproc++; if (iproc == nprocs) iproc = 0; - if (list[iproc] > 0) { + if (iproc == me) num_self = list[iproc]; + else if (list[iproc] > 0) { proc_send[isend] = iproc; num_send[isend] = list[iproc]; + list[iproc] = isend; isend++; } } + list[me] = 0; - // count = offsets into n-length index_send for each proc I send to + // count = offsets into index_send for each proc I send to + // m = ptr into index_self // index_send = list of which datums to send to each proc // 1st N1 values are datum indices for 1st proc, // next N2 values are datum indices for 2nd proc, etc - // offset_send = where each datum starts in send buffer count[0] = 0; for (i = 1; i < nsend; i++) count[i] = count[i-1] + num_send[i-1]; + m = 0; for (i = 0; i < n; i++) { - isend = list[proclist[i]]; - index_send[count[isend]++] = i; + iproc = proclist[i]; + if (iproc == me) index_self[m++] = i; + else { + isend = list[iproc]; + index_send[count[isend]++] = i; + } } // tell receivers how much data I send @@ -392,15 +487,16 @@ Irregular::PlanData *Irregular::create_data(int n, int *proclist, // num_recv = total size of message each proc sends me // nrecvsize = total size of data I recv - *nrecvsize = 0; + int nrecvsize = 0; for (i = 0; i < nrecv; i++) { MPI_Recv(&num_recv[i],1,MPI_INT,MPI_ANY_SOURCE,0,world,status); proc_recv[i] = status->MPI_SOURCE; - *nrecvsize += num_recv[i]; + nrecvsize += num_recv[i]; } + nrecvsize += num_self; // barrier to insure all MPI_ANY_SOURCE messages are received - // else another proc could proceed to irregular_perform() and send to me + // else another proc could proceed to exchange_data() and send to me MPI_Barrier(world); @@ -411,166 +507,105 @@ Irregular::PlanData *Irregular::create_data(int n, int *proclist, // initialize plan and return it - plan->nsend = nsend; - plan->nrecv = nrecv; - plan->sendmax = sendmax; + dplan->nsend = nsend; + dplan->nrecv = nrecv; + dplan->sendmax = sendmax; - plan->proc_send = proc_send; - plan->num_send = num_send; - plan->index_send = index_send; - plan->proc_recv = proc_recv; - plan->num_recv = num_recv; - plan->request = request; - plan->status = status; + dplan->proc_send = proc_send; + dplan->num_send = num_send; + dplan->index_send = index_send; + dplan->proc_recv = proc_recv; + dplan->num_recv = num_recv; + dplan->num_self = num_self; + dplan->index_self = index_self; - return plan; + dplan->request = request; + dplan->status = status; + + return nrecvsize; } /* ---------------------------------------------------------------------- - perform irregular communication of atoms - plan = previouly computed PlanAtom via create_atom() - sendbuf = list of atoms to send - sizes = # of doubles for each atom - recvbuf = received atoms -------------------------------------------------------------------------- */ - -void Irregular::exchange_atom(PlanAtom *plan, double *sendbuf, int *sizes, - double *recvbuf) -{ - int i,m,n,offset; - - // post all receives - - offset = 0; - for (int irecv = 0; irecv < plan->nrecv; irecv++) { - MPI_Irecv(&recvbuf[offset],plan->length_recv[irecv],MPI_DOUBLE, - plan->proc_recv[irecv],0,world,&plan->request[irecv]); - offset += plan->length_recv[irecv]; - } - - // allocate buf for largest send - - double *buf = (double *) memory->smalloc(plan->sendmax*sizeof(double), - "irregular:buf"); - - // send each message - // pack buf with list of datums (datum = one atom) - // m = index of datum in sendbuf - - n = 0; - for (int isend = 0; isend < plan->nsend; isend++) { - offset = 0; - for (i = 0; i < plan->num_send[isend]; i++) { - m = plan->index_send[n++]; - memcpy(&buf[offset],&sendbuf[plan->offset_send[m]], - sizes[m]*sizeof(double)); - offset += sizes[m]; - } - MPI_Send(buf,plan->length_send[isend],MPI_DOUBLE, - plan->proc_send[isend],0,world); - } - - // free temporary send buffer - - memory->sfree(buf); - - // wait on all incoming messages - - if (plan->nrecv) MPI_Waitall(plan->nrecv,plan->request,plan->status); -} - -/* ---------------------------------------------------------------------- - perform irregular communication of datums - plan = previouly computed PlanData via create_data() + communicate datums via PlanData sendbuf = list of datums to send nbytes = size of each datum - recvbuf = received datums + recvbuf = received datums (including copied from me) ------------------------------------------------------------------------- */ -void Irregular::exchange_data(PlanData *plan, char *sendbuf, int nbytes, - char *recvbuf) +void Irregular::exchange_data(char *sendbuf, int nbytes, char *recvbuf) { - int i,m,n,offset; + int i,m,n,offset,num_send; - // post all receives + // post all receives, starting after self copies - offset = 0; - for (int irecv = 0; irecv < plan->nrecv; irecv++) { - MPI_Irecv(&recvbuf[offset],nbytes*plan->num_recv[irecv],MPI_CHAR, - plan->proc_recv[irecv],0,world,&plan->request[irecv]); - offset += nbytes*plan->num_recv[irecv]; + offset = dplan->num_self*nbytes; + for (int irecv = 0; irecv < dplan->nrecv; irecv++) { + MPI_Irecv(&recvbuf[offset],dplan->num_recv[irecv]*nbytes,MPI_CHAR, + dplan->proc_recv[irecv],0,world,&dplan->request[irecv]); + offset += dplan->num_recv[irecv]*nbytes; } // allocate buf for largest send - char *buf = (char *) memory->smalloc(plan->sendmax*nbytes,"irregular:buf"); + char *buf = (char *) memory->smalloc(dplan->sendmax*nbytes,"irregular:buf"); // send each message - // pack buf with list of datums (datum = one atom) + // pack buf with list of datums // m = index of datum in sendbuf + int *index_send = dplan->index_send; + int nsend = dplan->nsend; n = 0; - for (int isend = 0; isend < plan->nsend; isend++) { - for (i = 0; i < plan->num_send[isend]; i++) { - m = plan->index_send[n++]; + + for (int isend = 0; isend < nsend; isend++) { + num_send = dplan->num_send[isend]; + for (i = 0; i < num_send; i++) { + m = index_send[n++]; memcpy(&buf[i*nbytes],&sendbuf[m*nbytes],nbytes); } - MPI_Send(buf,nbytes*plan->num_send[isend],MPI_CHAR, - plan->proc_send[isend],0,world); + MPI_Send(buf,dplan->num_send[isend]*nbytes,MPI_CHAR, + dplan->proc_send[isend],0,world); } // free temporary send buffer memory->sfree(buf); - // copy datums owned by self + // copy datums to self, put at beginning of recvbuf - for (i = 0; i < plan->num_self; i++) { - m = plan->index_self[i++]; + int *index_self = dplan->index_self; + int num_self = dplan->num_self; + + for (i = 0; i < num_self; i++) { + m = index_self[i]; memcpy(&recvbuf[i*nbytes],&sendbuf[m*nbytes],nbytes); } // wait on all incoming messages - if (plan->nrecv) MPI_Waitall(plan->nrecv,plan->request,plan->status); + if (dplan->nrecv) MPI_Waitall(dplan->nrecv,dplan->request,dplan->status); } /* ---------------------------------------------------------------------- - destroy an irregular communication plan for atoms + destroy communication plan for datums ------------------------------------------------------------------------- */ -void Irregular::destroy_atom(PlanAtom *plan) +void Irregular::destroy_data() { - delete [] plan->proc_send; - delete [] plan->length_send; - delete [] plan->num_send; - delete [] plan->index_send; - delete [] plan->offset_send; - delete [] plan->proc_recv; - delete [] plan->length_recv; - delete [] plan->request; - delete [] plan->status; - memory->sfree(plan); + delete [] dplan->proc_send; + delete [] dplan->num_send; + delete [] dplan->index_send; + delete [] dplan->proc_recv; + delete [] dplan->num_recv; + delete [] dplan->index_self; + delete [] dplan->request; + delete [] dplan->status; + memory->sfree(dplan); + dplan = NULL; } /* ---------------------------------------------------------------------- - destroy an irregular communication plan for datums -------------------------------------------------------------------------- */ - -void Irregular::destroy_data(PlanData *plan) -{ - delete [] plan->proc_send; - delete [] plan->num_send; - delete [] plan->index_send; - delete [] plan->proc_recv; - delete [] plan->num_recv; - delete [] plan->request; - delete [] plan->status; - memory->sfree(plan); -} - -/* ---------------------------------------------------------------------- - determine which proc owns atom with x coord + determine which proc owns atom with coord x[3] x will be in box (orthogonal) or lamda coords (triclinic) ------------------------------------------------------------------------- */ @@ -634,3 +669,13 @@ void Irregular::grow_recv(int n) "comm:buf_recv"); } +/* ---------------------------------------------------------------------- + return # of bytes of allocated memory +------------------------------------------------------------------------- */ + +double Irregular::memory_usage() +{ + double bytes = maxsend * sizeof(double); + bytes += maxrecv * sizeof(double); + return bytes; +} diff --git a/src/irregular.h b/src/irregular.h index 7306a58606..f535e10527 100644 --- a/src/irregular.h +++ b/src/irregular.h @@ -20,21 +20,45 @@ namespace LAMMPS_NS { class Irregular : protected Pointers { public: - struct PlanAtom { // plan for irregular communication of atoms + Irregular(class LAMMPS *); + ~Irregular(); + void migrate_atoms(); + int create_data(int, int *); + void exchange_data(char *, int, char *); + void destroy_data(); + double memory_usage(); + + private: + int me,nprocs; + int triclinic; + int map_style; + int *procgrid; + int ***grid2proc; + + int maxsend,maxrecv; // size of buffers in # of doubles + double *buf_send,*buf_recv; + + // plan for irregular communication of atoms + // no params refer to atoms copied to self + + struct PlanAtom { int nsend; // # of messages to send int nrecv; // # of messages to recv int sendmax; // # of doubles in largest send message int *proc_send; // procs to send to int *length_send; // # of doubles to send to each proc - int *num_send; // # of datums to send to each proc - int *index_send; // list of which datums to send to each proc - int *offset_send; // where each datum starts in send buffer + int *num_send; // # of atoms to send to each proc + int *index_send; // list of which atoms to send to each proc + int *offset_send; // where each atom starts in send buffer int *proc_recv; // procs to recv from int *length_recv; // # of doubles to recv from each proc MPI_Request *request; // MPI requests for posted recvs MPI_Status *status; // MPI statuses for WaitAll }; + // plan for irregular communication of datums + // only 2 self params refer to atoms copied to self + struct PlanData { // plan for irregular communication of data int nsend; // # of messages to send int nrecv; // # of messages to recv @@ -43,32 +67,19 @@ class Irregular : protected Pointers { int *num_send; // # of datums to send to each proc int *index_send; // list of which datums to send to each proc int *proc_recv; // procs to recv from - int *num_recv; // # of doubles to recv from each proc - int num_self; - int *index_self; + int *num_recv; // # of datums to recv from each proc + int num_self; // # of datums to copy to self + int *index_self; // list of which datums to copy to self MPI_Request *request; // MPI requests for posted recvs MPI_Status *status; // MPI statuses for WaitAll }; - Irregular(class LAMMPS *); - ~Irregular(); - void migrate_atoms(); - struct PlanData *create_data(int, int *, int *); - void exchange_data(PlanData *, char *, int, char *); - void destroy_data(PlanData *); + PlanAtom *aplan; + PlanData *dplan; - private: - int me,nprocs; - int triclinic; - int map_style; - int *procgrid; - int ***grid2proc; - int maxsend,maxrecv; - double *buf_send,*buf_recv; - - struct PlanAtom *create_atom(int, int *, int *, int *); - void exchange_atom(PlanAtom *, double *, int *, double *); - void destroy_atom(PlanAtom *); + int create_atom(int, int *, int *); + void exchange_atom(double *, int *, double *); + void destroy_atom(); int coord2proc(double *); void grow_send(int,int); // reallocate send buffer