diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp index ccdcf9883c..3127621fef 100644 --- a/src/GRANULAR/pair_gran_hooke_history.cpp +++ b/src/GRANULAR/pair_gran_hooke_history.cpp @@ -466,6 +466,11 @@ double PairGranHookeHistory::init_one(int i, int j) void PairGranHookeHistory::write_restart(FILE *fp) { write_restart_settings(fp); + + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) + fwrite(&setflag[i][j],sizeof(int),1,fp); } /* ---------------------------------------------------------------------- @@ -476,6 +481,14 @@ void PairGranHookeHistory::read_restart(FILE *fp) { read_restart_settings(fp); allocate(); + + int i,j; + int me = comm->me; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); + MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + } } /* ---------------------------------------------------------------------- diff --git a/src/fix_read_restart.cpp b/src/fix_read_restart.cpp new file mode 100644 index 0000000000..3ed607d683 --- /dev/null +++ b/src/fix_read_restart.cpp @@ -0,0 +1,126 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "stdlib.h" +#include "fix_read_restart.h" +#include "atom.h" +#include "memory.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +FixReadRestart::FixReadRestart(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + nextra = atoi(arg[3]); + int nfix = atoi(arg[4]); + + // perform initial allocation of atom-based array + // register with Atom class + + count = NULL; + extra = NULL; + grow_arrays(atom->nmax); + atom->add_callback(0); + + // extra = copy of atom->extra + + double **atom_extra = atom->extra; + int nlocal = atom->nlocal; + int i,j,m; + + for (i = 0; i < nlocal; i++) { + m = 0; + for (j = 0; j < nfix; j++) m += static_cast (atom_extra[i][m]); + count[i] = m; + for (j = 0; j < m; j++) extra[i][j] = atom_extra[i][j]; + } +} + +/* ---------------------------------------------------------------------- */ + +FixReadRestart::~FixReadRestart() +{ + // unregister callback to this fix from Atom class + + atom->delete_callback(id,0); + + // delete locally stored arrays + + memory->sfree(count); + memory->destroy_2d_double_array(extra); +} + +/* ---------------------------------------------------------------------- */ + +int FixReadRestart::setmask() +{ + int mask = 0; + return mask; +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based array +------------------------------------------------------------------------- */ + +double FixReadRestart::memory_usage() +{ + double bytes = atom->nmax*nextra * sizeof(double); + bytes += atom->nmax * sizeof(int); + return bytes; +} + +/* ---------------------------------------------------------------------- + allocate atom-based array +------------------------------------------------------------------------- */ + +void FixReadRestart::grow_arrays(int nmax) +{ + count = + (int *) memory->srealloc(count,nmax*sizeof(int),"read_restart:count"); + extra = + memory->grow_2d_double_array(extra,nmax,nextra,"read_restart:extra"); +} + +/* ---------------------------------------------------------------------- + copy values within local atom-based array +------------------------------------------------------------------------- */ + +void FixReadRestart::copy_arrays(int i, int j) +{ + count[j] = count[i]; + for (int m = 0; m < count[i]; m++) extra[j][m] = extra[i][m]; +} + +/* ---------------------------------------------------------------------- + pack values in local atom-based array for exchange with another proc +------------------------------------------------------------------------- */ + +int FixReadRestart::pack_exchange(int i, double *buf) +{ + buf[0] = count[i]; + for (int m = 0; m < count[i]; m++) buf[m+1] = extra[i][m]; + return count[i]+1; +} + +/* ---------------------------------------------------------------------- + unpack values in local atom-based array from exchange with another proc +------------------------------------------------------------------------- */ + +int FixReadRestart::unpack_exchange(int nlocal, double *buf) +{ + count[nlocal] = static_cast (buf[0]); + for (int m = 0; m < count[nlocal]; m++) extra[nlocal][m] = buf[m+1]; + return count[nlocal]+1; +} diff --git a/src/fix_read_restart.h b/src/fix_read_restart.h new file mode 100644 index 0000000000..d1b985b4bb --- /dev/null +++ b/src/fix_read_restart.h @@ -0,0 +1,49 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(READ_RESTART,FixReadRestart) + +#else + +#ifndef LMP_FIX_READ_RESTART_H +#define LMP_FIX_READ_RESTART_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixReadRestart : public Fix { + public: + int *count; + double **extra; + + FixReadRestart(class LAMMPS *, int, char **); + ~FixReadRestart(); + int setmask(); + + double memory_usage(); + void grow_arrays(int); + void copy_arrays(int, int); + int pack_exchange(int, double *); + int unpack_exchange(int, double *); + + private: + int nextra; // max number of extra values for any atom +}; + +} + +#endif +#endif diff --git a/src/modify.cpp b/src/modify.cpp index 2b39865854..139a87161e 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -91,6 +91,7 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp) nfix_restart_peratom = 0; id_restart_peratom = style_restart_peratom = NULL; index_restart_peratom = NULL; + nfix_restart_save = 0; ncompute = maxcompute = 0; compute = NULL; @@ -146,8 +147,9 @@ void Modify::init() int i,j; // delete storage of restart info since it is not valid after 1st run + // read_restart sets nfix_restart_save so will persist to actual 1st run - restart_deallocate(); + if (!nfix_restart_save) restart_deallocate(); // create lists of fixes to call at each stage of run diff --git a/src/modify.h b/src/modify.h index 5fb936a86b..0cdf25600f 100644 --- a/src/modify.h +++ b/src/modify.h @@ -31,7 +31,8 @@ class Modify : protected Pointers { int restart_pbc_any; // 1 if any fix sets restart_pbc int nfix_restart_global; // stored fix global info from restart file - int nfix_restart_peratom; // stored fix peratom info from restart file + int nfix_restart_peratom; // stored fix peratom info from restart file + int nfix_restart_save; // 1 if init() should not whack restart info class Fix **fix; // list of fixes int *fmask; // bit mask for when each fix is applied diff --git a/src/read_restart.cpp b/src/read_restart.cpp index bf908caaba..c9356d0edb 100644 --- a/src/read_restart.cpp +++ b/src/read_restart.cpp @@ -23,6 +23,8 @@ #include "comm.h" #include "update.h" #include "modify.h" +#include "fix.h" +#include "fix_read_restart.h" #include "group.h" #include "force.h" #include "pair.h" @@ -70,6 +72,7 @@ void ReadRestart::command(int narg, char **arg) error->all("Cannot read_restart after simulation box is defined"); MPI_Comm_rank(world,&me); + MPI_Comm_size(world,&nprocs); // if filename contains "*", search dir for latest restart file @@ -112,8 +115,8 @@ void ReadRestart::command(int narg, char **arg) // problem setup using info from header int n; - if (comm->nprocs == 1) n = static_cast (atom->natoms); - else n = static_cast (LB_FACTOR * atom->natoms / comm->nprocs); + if (nprocs == 1) n = static_cast (atom->natoms); + else n = static_cast (LB_FACTOR * atom->natoms / nprocs); atom->allocate_type_arrays(); atom->avec->grow(n); @@ -136,85 +139,161 @@ void ReadRestart::command(int narg, char **arg) atom->nextra_store = nextra; atom->extra = memory->create_2d_double_array(n,nextra,"atom:extra"); - // if single file: - // proc 0 reads atoms from file, one chunk per proc (nprocs_file) - // else if one file per proc: - // proc 0 reads chunks from series of files (nprocs_file) - // proc 0 bcasts each chunk to other procs + // single file: + // nprocs_file = # of chunks in file + // proc 0 reads chunks one at a time and bcasts it to other procs // each proc unpacks the atoms, saving ones in it's sub-domain - // check for atom in sub-domain is different for orthogonal vs triclinic box + // check for atom in sub-domain differs for orthogonal vs triclinic box + // close restart file when done AtomVec *avec = atom->avec; - int triclinic = domain->triclinic; int maxbuf = 0; double *buf = NULL; - - double *x,lamda[3]; - double *coord,*sublo,*subhi; - if (triclinic == 0) { - sublo = domain->sublo; - subhi = domain->subhi; - } else { - sublo = domain->sublo_lamda; - subhi = domain->subhi_lamda; - } - int m; - char *perproc = new char[strlen(file) + 16]; - char *ptr = strchr(file,'%'); - for (int iproc = 0; iproc < nprocs_file; iproc++) { - if (me == 0) { - if (multiproc) { - fclose(fp); - *ptr = '\0'; - sprintf(perproc,"%s%d%s",file,iproc,ptr+1); - *ptr = '%'; - fp = fopen(perproc,"rb"); - if (fp == NULL) { - char str[128]; - sprintf(str,"Cannot open restart file %s",perproc); - error->one(str); - } + if (multiproc == 0) { + int triclinic = domain->triclinic; + double *x,lamda[3]; + double *coord,*sublo,*subhi; + if (triclinic == 0) { + sublo = domain->sublo; + subhi = domain->subhi; + } else { + sublo = domain->sublo_lamda; + subhi = domain->subhi_lamda; + } + + for (int iproc = 0; iproc < nprocs_file; iproc++) { + if (me == 0) fread(&n,sizeof(int),1,fp); + MPI_Bcast(&n,1,MPI_INT,0,world); + if (n > maxbuf) { + maxbuf = n; + memory->sfree(buf); + buf = (double *) memory->smalloc(maxbuf*sizeof(double), + "read_restart:buf"); } + + if (n > 0) { + if (me == 0) fread(buf,sizeof(double),n,fp); + MPI_Bcast(buf,n,MPI_DOUBLE,0,world); + } + + m = 0; + while (m < n) { + x = &buf[m+1]; + if (triclinic) { + domain->x2lamda(x,lamda); + coord = lamda; + } else coord = x; + + if (coord[0] >= sublo[0] && coord[0] < subhi[0] && + coord[1] >= sublo[1] && coord[1] < subhi[1] && + coord[2] >= sublo[2] && coord[2] < subhi[2]) { + m += avec->unpack_restart(&buf[m]); + } + else m += static_cast (buf[m]); + } + } + + if (me == 0) fclose(fp); + + // one file per proc: + // nprocs_file = # of files + // each proc reads 1/P fraction of files, keeping all atoms in the files + // perform irregular comm to migrate atoms to correct procs + // close restart file when done + + } else { + if (me == 0) fclose(fp); + char *perproc = new char[strlen(file) + 16]; + char *ptr = strchr(file,'%'); + + for (int iproc = me; iproc < nprocs_file; iproc += nprocs) { + *ptr = '\0'; + sprintf(perproc,"%s%d%s",file,iproc,ptr+1); + *ptr = '%'; + fp = fopen(perproc,"rb"); + if (fp == NULL) { + char str[128]; + sprintf(str,"Cannot open restart file %s",perproc); + error->one(str); + } + fread(&n,sizeof(int),1,fp); + if (n > maxbuf) { + maxbuf = n; + memory->sfree(buf); + buf = (double *) memory->smalloc(maxbuf*sizeof(double), + "read_restart:buf"); + } + if (n > 0) fread(buf,sizeof(double),n,fp); + + m = 0; + while (m < n) m += avec->unpack_restart(&buf[m]); + fclose(fp); } - MPI_Bcast(&n,1,MPI_INT,0,world); - if (n > maxbuf) { - maxbuf = n; - delete [] buf; - buf = new double[maxbuf]; + delete [] perproc; + + // save modify->nfix_restart values, so Modify::init() won't whack them + // create a temporary fix to hold extra atom info + // necessary b/c Atom::init() will destroy atom->extra + + if (nextra) { + modify->nfix_restart_save = 1; + char cextra[8],fixextra[8]; + sprintf(cextra,"%d",nextra); + sprintf(fixextra,"%d",modify->nfix_restart_peratom); + char **newarg = new char*[5]; + newarg[0] = (char *) "_read_restart"; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "READ_RESTART"; + newarg[3] = cextra; + newarg[4] = fixextra; + modify->add_fix(5,newarg); + delete [] newarg; } - if (n > 0) { - if (me == 0) fread(buf,sizeof(double),n,fp); - MPI_Bcast(buf,n,MPI_DOUBLE,0,world); - } + // init entire system before comm->irregular + // comm::init needs neighbor::init needs pair::init needs kspace::init, etc + // move atoms to new processors via irregular() + // in case read by totally different proc than wrote restart file - m = 0; - while (m < n) { - x = &buf[m+1]; - if (triclinic) { - domain->x2lamda(x,lamda); - coord = lamda; - } else coord = x; + if (me == 0 && screen) + fprintf(screen," system init for parallel read_restart ...\n"); + lmp->init(); + if (domain->triclinic) domain->x2lamda(atom->nlocal); + domain->reset_box(); + comm->setup(); + comm->irregular(); + if (domain->triclinic) domain->lamda2x(atom->nlocal); - if (coord[0] >= sublo[0] && coord[0] < subhi[0] && - coord[1] >= sublo[1] && coord[1] < subhi[1] && - coord[2] >= sublo[2] && coord[2] < subhi[2]) - m += avec->unpack_restart(&buf[m]); - else m += static_cast (buf[m]); + // unsave modify->nfix_restart values and atom->extra + // destroy temporary fix + + if (nextra) { + modify->nfix_restart_save = 0; + atom->nextra_store = nextra; + atom->extra = memory->create_2d_double_array(atom->nmax,nextra, + "atom:extra"); + int ifix = modify->find_fix("_read_restart"); + FixReadRestart *fix = (FixReadRestart *) modify->fix[ifix]; + int *count = fix->count; + double **extra = fix->extra; + double **atom_extra = atom->extra; + int nlocal = atom->nlocal; + for (int i = 0; i < nlocal; i++) + for (int j = 0; j < count[i]; j++) + atom_extra[i][j] = extra[i][j]; + modify->delete_fix("_read_restart"); } } - // close restart file and clean-up memory - - if (me == 0) fclose(fp); - delete [] buf; + // clean-up memory + delete [] file; - delete [] perproc; + memory->sfree(buf); // check that all atoms were assigned to procs @@ -422,9 +501,9 @@ void ReadRestart::header() pz != comm->user_procgrid[2]) && me == 0) error->warning("Restart file used different 3d processor grid"); - // don't set newton_pair, leave input script value unchanged - // set newton_bond from restart file - // warn if different and input script settings are not default + // don't set newton_pair, leave input script value unchanged + // set newton_bond from restart file + // warn if different and input script settings are not default } else if (flag == NEWTON_PAIR) { int newton_pair_file = read_int(); diff --git a/src/read_restart.h b/src/read_restart.h index 2db6421d95..fd35f8c97b 100644 --- a/src/read_restart.h +++ b/src/read_restart.h @@ -31,9 +31,9 @@ class ReadRestart : protected Pointers { void command(int, char **); private: - int me; + int me,nprocs,nprocs_file; FILE *fp; - int nprocs_file; + int nfix_restart_global,nfix_restart_peratom; void file_search(char *, char *); void header();