git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@4774 f3b2605a-c512-4ea7-a41b-209d697bcdaa

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);
+    }
 }
 
 /* ----------------------------------------------------------------------

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<int> (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<int> (buf[0]);
+  for (int m = 0; m < count[nlocal]; m++) extra[nlocal][m] = buf[m+1];
+  return count[nlocal]+1;
+}

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

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
 

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

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<int> (atom->natoms);
+  if (nprocs == 1) n = static_cast<int> (atom->natoms);
-  else n = static_cast<int> (LB_FACTOR * atom->natoms / comm->nprocs);
+  else n = static_cast<int> (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:
+  // single file:
-  //   proc 0 reads atoms from file, one chunk per proc (nprocs_file)
+  // nprocs_file = # of chunks in file
-  // else if one file per proc:
+  // proc 0 reads chunks one at a time and bcasts it to other procs
-  //   proc 0 reads chunks from series of files (nprocs_file)
-  // proc 0 bcasts each chunk 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 (multiproc == 0) {
-    if (me == 0) {
+    int triclinic = domain->triclinic;
-      if (multiproc) {
+    double *x,lamda[3];
-	fclose(fp);
+    double *coord,*sublo,*subhi;
-	*ptr = '\0';
+    if (triclinic == 0) {
-	sprintf(perproc,"%s%d%s",file,iproc,ptr+1);
+      sublo = domain->sublo;
-	*ptr = '%';
+      subhi = domain->subhi;
-	fp = fopen(perproc,"rb");
+    } else {
-	if (fp == NULL) {
+      sublo = domain->sublo_lamda;
-	  char str[128];
+      subhi = domain->subhi_lamda;
-	  sprintf(str,"Cannot open restart file %s",perproc);
+    }
-	  error->one(str);
+
-	}
+    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<int> (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);
+    delete [] perproc;
-    if (n > maxbuf) {
+
-      maxbuf = n;
+    // save modify->nfix_restart values, so Modify::init() won't whack them
-      delete [] buf;
+    // create a temporary fix to hold extra atom info
-      buf = new double[maxbuf];
+    // 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) {
+    // init entire system before comm->irregular
-      if (me == 0) fread(buf,sizeof(double),n,fp);
+    // comm::init needs neighbor::init needs pair::init needs kspace::init, etc
-      MPI_Bcast(buf,n,MPI_DOUBLE,0,world);
+    // move atoms to new processors via irregular()
-    }
+    // in case read by totally different proc than wrote restart file
 
-    m = 0;
+    if (me == 0 && screen)
-    while (m < n) {
+      fprintf(screen,"  system init for parallel read_restart ...\n");
-      x = &buf[m+1];
+    lmp->init();
-      if (triclinic) {
+    if (domain->triclinic) domain->x2lamda(atom->nlocal);
-	domain->x2lamda(x,lamda);
+    domain->reset_box();
-	coord = lamda;
+    comm->setup();
-      } else coord = x;
+    comm->irregular();
+    if (domain->triclinic) domain->lamda2x(atom->nlocal);
 
-      if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
+    // unsave modify->nfix_restart values and atom->extra
-	  coord[1] >= sublo[1] && coord[1] < subhi[1] &&
+    // destroy temporary fix
-	  coord[2] >= sublo[2] && coord[2] < subhi[2])
+
-	m += avec->unpack_restart(&buf[m]);
+    if (nextra) {
-      else m += static_cast<int> (buf[m]);
+      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
+  // clean-up memory
-  
+
-  if (me == 0) fclose(fp);
-  delete [] buf;
   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
+    // don't set newton_pair, leave input script value unchanged
-      // set newton_bond from restart file
+    // set newton_bond from restart file
-      // warn if different and input script settings are not default
+    // warn if different and input script settings are not default
 
     } else if (flag == NEWTON_PAIR) {
       int newton_pair_file = read_int();

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();