diff --git a/src/compute_property_local.cpp b/src/compute_property_local.cpp index 22a2228292..61646ee0f7 100644 --- a/src/compute_property_local.cpp +++ b/src/compute_property_local.cpp @@ -20,6 +20,9 @@ #include "memory.h" #include "error.h" + +#include "comm.h" + using namespace LAMMPS_NS; // customize by adding keyword @@ -174,6 +177,7 @@ void ComputePropertyLocal::init() else if (kindflag == IMPROPER) ncount = count_impropers(0); if (ncount > nmax) reallocate(ncount); + size_local_rows = ncount; } /* ---------------------------------------------------------------------- */ @@ -190,6 +194,7 @@ void ComputePropertyLocal::compute_local() else if (kindflag == IMPROPER) ncount = count_impropers(0); if (ncount > nmax) reallocate(ncount); + size_local_rows = ncount; if (kindflag == BOND) ncount = count_bonds(1); else if (kindflag == ANGLE) ncount = count_angles(1); @@ -245,6 +250,7 @@ int ComputePropertyLocal::count_bonds(int flag) m++; } } + return m; } @@ -285,6 +291,7 @@ int ComputePropertyLocal::count_angles(int flag) m++; } } + return m; } @@ -328,6 +335,7 @@ int ComputePropertyLocal::count_dihedrals(int flag) m++; } } + return m; } @@ -371,6 +379,7 @@ int ComputePropertyLocal::count_impropers(int flag) m++; } } + return m; } @@ -380,8 +389,6 @@ void ComputePropertyLocal::reallocate(int n) { // grow vector or array and indices array - size_local_rows = n; - while (nmax < n) nmax += DELTA; if (nvalues == 1) { memory->sfree(vector); diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp index ef455442e6..0532f988e1 100644 --- a/src/dump_custom.cpp +++ b/src/dump_custom.cpp @@ -65,8 +65,8 @@ DumpCustom::DumpCustom(LAMMPS *lmp, int narg, char **arg) : // computes, fixes, variables which the dump accesses - field2index = (int *) memory->smalloc(nfield*sizeof(int),"dump:field2index"); - argindex = (int *) memory->smalloc(nfield*sizeof(int),"dump:argindex"); + field2index = new int[nfield]; + argindex = new int[nfield]; ncompute = 0; id_compute = NULL; @@ -122,14 +122,13 @@ DumpCustom::~DumpCustom() { delete [] pack_choice; delete [] vtype; + delete [] field2index; + delete [] argindex; memory->sfree(thresh_array); memory->sfree(thresh_op); memory->sfree(thresh_value); - memory->sfree(field2index); - memory->sfree(argindex); - for (int i = 0; i < ncompute; i++) delete [] id_compute[i]; memory->sfree(id_compute); delete [] compute; diff --git a/src/dump_local.cpp b/src/dump_local.cpp new file mode 100644 index 0000000000..723bafb49a --- /dev/null +++ b/src/dump_local.cpp @@ -0,0 +1,427 @@ +/* ---------------------------------------------------------------------- + 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 "mpi.h" +#include "string.h" +#include "stdlib.h" +#include "dump_local.h" +#include "atom.h" +#include "modify.h" +#include "fix.h" +#include "compute.h" +#include "domain.h" +#include "update.h" +#include "memory.h" +#include "error.h" + + +#include "comm.h" + +using namespace LAMMPS_NS; + +enum{INDEX,COMPUTE,FIX}; +enum{INT,DOUBLE}; + +#define INVOKED_LOCAL 16 + +/* ---------------------------------------------------------------------- */ + +DumpLocal::DumpLocal(LAMMPS *lmp, int narg, char **arg) : + Dump(lmp, narg, arg) +{ + if (narg == 5) error->all("No dump local arguments specified"); + + nevery = atoi(arg[3]); + + size_one = nfield = narg-5; + pack_choice = new FnPtrPack[nfield]; + vtype = new int[nfield]; + + // computes & fixes which the dump accesses + + field2index = new int[nfield]; + argindex = new int[nfield]; + + ncompute = 0; + id_compute = NULL; + compute = NULL; + + nfix = 0; + id_fix = NULL; + fix = NULL; + + // process keywords + + parse_fields(narg,arg); + + // setup format strings + + vformat = new char*[size_one]; + + format_default = new char[3*size_one+1]; + format_default[0] = '\0'; + + for (int i = 0; i < size_one; i++) { + if (vtype[i] == INT) format_default = strcat(format_default,"%d "); + else format_default = strcat(format_default,"%g "); + vformat[i] = NULL; + } + + // setup column string + + int n = 0; + for (int iarg = 5; iarg < narg; iarg++) n += strlen(arg[iarg]) + 2; + columns = new char[n]; + columns[0] = '\0'; + for (int iarg = 5; iarg < narg; iarg++) { + strcat(columns,arg[iarg]); + strcat(columns," "); + } + + label = "BONDS"; +} + +/* ---------------------------------------------------------------------- */ + +DumpLocal::~DumpLocal() +{ + delete [] pack_choice; + delete [] vtype; + delete [] field2index; + delete [] argindex; + + for (int i = 0; i < ncompute; i++) delete [] id_compute[i]; + memory->sfree(id_compute); + delete [] compute; + + for (int i = 0; i < nfix; i++) delete [] id_fix[i]; + memory->sfree(id_fix); + delete [] fix; + + for (int i = 0; i < size_one; i++) delete [] vformat[i]; + delete [] vformat; + + delete [] columns; +} + +/* ---------------------------------------------------------------------- */ + +void DumpLocal::init() +{ + delete [] format; + char *str; + if (format_user) str = format_user; + else str = format_default; + + int n = strlen(str) + 1; + format = new char[n]; + strcpy(format,str); + + // tokenize the format string and add space at end of each format element + + char *ptr; + for (int i = 0; i < size_one; i++) { + if (i == 0) ptr = strtok(format," \0"); + else ptr = strtok(NULL," \0"); + delete [] vformat[i]; + vformat[i] = new char[strlen(ptr) + 2]; + strcpy(vformat[i],ptr); + vformat[i] = strcat(vformat[i]," "); + } + + // find current ptr for each compute,fix,variable + // check that fix frequency is acceptable + + int icompute; + for (int i = 0; i < ncompute; i++) { + icompute = modify->find_compute(id_compute[i]); + if (icompute < 0) error->all("Could not find dump custom compute ID"); + compute[i] = modify->compute[icompute]; + } + + int ifix; + for (int i = 0; i < nfix; i++) { + ifix = modify->find_fix(id_fix[i]); + if (ifix < 0) error->all("Could not find dump custom fix ID"); + fix[i] = modify->fix[ifix]; + if (nevery % modify->fix[ifix]->peratom_freq) + error->all("Dump custom and fix not computed at compatible times"); + } + + // open single file, one time only + + if (multifile == 0) openfile(); +} + +/* ---------------------------------------------------------------------- */ + +void DumpLocal::write_header(int ndump) +{ + if (me == 0) { + fprintf(fp,"ITEM: TIMESTEP\n"); + fprintf(fp,"%d\n",update->ntimestep); + fprintf(fp,"ITEM: NUMBER OF %s\n",label); + fprintf(fp,"%d\n",ndump); + fprintf(fp,"ITEM: %s\n",label); + } +} + +/* ---------------------------------------------------------------------- */ + +int DumpLocal::count() +{ + int i; + + // invoke Computes for local quantities + + if (ncompute) { + int ntimestep = update->ntimestep; + for (i = 0; i < ncompute; i++) { + if (!(compute[i]->invoked_flag & INVOKED_LOCAL)) { + compute[i]->compute_local(); + compute[i]->invoked_flag |= INVOKED_LOCAL; + nmine = compute[i]->size_local_rows; + } + } + } + + return nmine; +} + +/* ---------------------------------------------------------------------- */ + +int DumpLocal::pack() +{ + for (int n = 0; n < size_one; n++) (this->*pack_choice[n])(n); + return nmine*size_one; +} + +/* ---------------------------------------------------------------------- */ + +void DumpLocal::write_data(int n, double *buf) +{ + int i,j; + + int m = 0; + for (i = 0; i < n; i++) { + for (j = 0; j < size_one; j++) { + if (vtype[j] == INT) fprintf(fp,vformat[j],static_cast (buf[m])); + else fprintf(fp,vformat[j],buf[m]); + m++; + } + fprintf(fp,"\n"); + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpLocal::parse_fields(int narg, char **arg) +{ + // customize by adding to if statement + + int i; + for (int iarg = 5; iarg < narg; iarg++) { + i = iarg-5; + + if (strcmp(arg[iarg],"index") == 0) { + pack_choice[i] = &DumpLocal::pack_index; + vtype[i] = INT; + + // compute value = c_ID + // if no trailing [], then arg is set to 0, else arg is int between [] + + } else if (strncmp(arg[iarg],"c_",2) == 0) { + pack_choice[i] = &DumpLocal::pack_compute; + vtype[i] = DOUBLE; + + int n = strlen(arg[iarg]); + char *suffix = new char[n]; + strcpy(suffix,&arg[iarg][2]); + + char *ptr = strchr(suffix,'['); + if (ptr) { + if (suffix[strlen(suffix)-1] != ']') + error->all("Invalid keyword in dump local command"); + argindex[i] = atoi(ptr+1); + *ptr = '\0'; + } else argindex[i] = 0; + + n = modify->find_compute(suffix); + if (n < 0) error->all("Could not find dump local compute ID"); + if (modify->compute[n]->local_flag == 0) + error->all("Dump local compute does not compute local info"); + if (argindex[i] == 0 && modify->compute[n]->size_local_cols > 0) + error->all("Dump local compute does not calculate local vector"); + if (argindex[i] > 0 && modify->compute[n]->size_local_cols == 0) + error->all("Dump local compute does not calculate local array"); + if (argindex[i] > 0 && + argindex[i] > modify->compute[n]->size_local_cols) + error->all("Dump local compute vector is accessed out-of-range"); + + field2index[i] = add_compute(suffix); + delete [] suffix; + + // fix value = f_ID + // if no trailing [], then arg is set to 0, else arg is between [] + + } else if (strncmp(arg[iarg],"f_",2) == 0) { + pack_choice[i] = &DumpLocal::pack_fix; + vtype[i] = DOUBLE; + + int n = strlen(arg[iarg]); + char *suffix = new char[n]; + strcpy(suffix,&arg[iarg][2]); + + char *ptr = strchr(suffix,'['); + if (ptr) { + if (suffix[strlen(suffix)-1] != ']') + error->all("Invalid keyword in dump local command"); + argindex[i] = atoi(ptr+1); + *ptr = '\0'; + } else argindex[i] = 0; + + n = modify->find_fix(suffix); + if (n < 0) error->all("Could not find dump local fix ID"); + if (modify->fix[n]->local_flag == 0) + error->all("Dump local fix does not compute local info"); + if (argindex[i] == 0 && modify->fix[n]->size_local_cols > 0) + error->all("Dump local fix does not compute local vector"); + if (argindex[i] > 0 && modify->fix[n]->size_local_cols == 0) + error->all("Dump local fix does not compute local array"); + if (argindex[i] > 0 && + argindex[i] > modify->fix[n]->size_local_cols) + error->all("Dump local fix vector is accessed out-of-range"); + + field2index[i] = add_fix(suffix); + delete [] suffix; + + } else error->all("Invalid keyword in dump local command"); + } +} + +/* ---------------------------------------------------------------------- + add Compute to list of Compute objects used by dump + return index of where this Compute is in list + if already in list, do not add, just return index, else add to list +------------------------------------------------------------------------- */ + +int DumpLocal::add_compute(char *id) +{ + int icompute; + for (icompute = 0; icompute < ncompute; icompute++) + if (strcmp(id,id_compute[icompute]) == 0) break; + if (icompute < ncompute) return icompute; + + id_compute = (char **) + memory->srealloc(id_compute,(ncompute+1)*sizeof(char *),"dump:id_compute"); + delete [] compute; + compute = new Compute*[ncompute+1]; + + int n = strlen(id) + 1; + id_compute[ncompute] = new char[n]; + strcpy(id_compute[ncompute],id); + ncompute++; + return ncompute-1; +} + +/* ---------------------------------------------------------------------- + add Fix to list of Fix objects used by dump + return index of where this Fix is in list + if already in list, do not add, just return index, else add to list +------------------------------------------------------------------------- */ + +int DumpLocal::add_fix(char *id) +{ + int ifix; + for (ifix = 0; ifix < nfix; ifix++) + if (strcmp(id,id_fix[ifix]) == 0) break; + if (ifix < nfix) return ifix; + + id_fix = (char **) + memory->srealloc(id_fix,(nfix+1)*sizeof(char *),"dump:id_fix"); + delete [] fix; + fix = new Fix*[nfix+1]; + + int n = strlen(id) + 1; + id_fix[nfix] = new char[n]; + strcpy(id_fix[nfix],id); + nfix++; + return nfix-1; +} + +/* ---------------------------------------------------------------------- + extraction of Compute, Fix results +------------------------------------------------------------------------- */ + +void DumpLocal::pack_compute(int n) +{ + double *vector = compute[field2index[n]]->vector_local; + double **array = compute[field2index[n]]->array_local; + int ncount = compute[field2index[n]]->size_local_rows; + int index = argindex[n]; + + if (index == 0) { + for (int i = 0; i < ncount; i++) { + buf[n] = vector[i]; + n += size_one; + } + } else { + index--; + for (int i = 0; i < ncount; i++) { + buf[n] = array[i][index]; + n += size_one; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpLocal::pack_fix(int n) +{ + double *vector = fix[field2index[n]]->vector_local; + double **array = fix[field2index[n]]->array_local; + int index = argindex[n]; + + if (index == 0) { + for (int i = 0; i < nmine; i++) { + buf[n] = vector[i]; + n += size_one; + } + } else { + index--; + for (int i = 0; i < nmine; i++) { + buf[n] = array[i][index]; + n += size_one; + } + } +} + +/* ---------------------------------------------------------------------- + one method for every keyword dump local can output + the local value is packed into buf starting at n with stride size_one + customize a new keyword by adding a method +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- */ + +void DumpLocal::pack_index(int n) +{ + int index; + MPI_Scan(&nmine,&index,1,MPI_INT,MPI_SUM,world); + index -= nmine; + + for (int i = 0; i < nmine; i++) { + buf[n] = ++index; + n += size_one; + } +} diff --git a/src/dump_local.h b/src/dump_local.h new file mode 100644 index 0000000000..a4fc4fe141 --- /dev/null +++ b/src/dump_local.h @@ -0,0 +1,74 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifndef DUMP_LOCAL_H +#define DUMP_LOCAL_H + +#include "dump.h" + +namespace LAMMPS_NS { + +class DumpLocal : public Dump { + public: + DumpLocal(LAMMPS *, int, char **); + ~DumpLocal(); + void init(); + + private: + int nevery; // dump frequency to check Fix against + char *label; // string for dump file header + + int nmine; // # of lines I am dumping + int *vtype; // type of each vector (INT, DOUBLE) + char **vformat; // format string for each vector element + + char *columns; // column labels + + int nfield; // # of keywords listed by user + + int *field2index; // which compute,fix,variable calcs this field + int *argindex; // index into compute,fix scalar_atom,vector_atom + // 0 for scalar_atom, 1-N for vector_atom values + + int ncompute; // # of Compute objects used by dump + char **id_compute; // their IDs + class Compute **compute; // list of ptrs to the Compute objects + + int nfix; // # of Fix objects used by dump + char **id_fix; // their IDs + class Fix **fix; // list of ptrs to the Fix objects + + // private methods + + void write_header(int); + int count(); + int pack(); + void write_data(int, double *); + + void parse_fields(int, char **); + int add_compute(char *); + int add_fix(char *); + + // customize by adding a method prototype + + typedef void (DumpLocal::*FnPtrPack)(int); + FnPtrPack *pack_choice; // ptrs to pack functions + + void pack_index(int); + void pack_compute(int); + void pack_fix(int); +}; + +} + +#endif diff --git a/src/output.cpp b/src/output.cpp index 49784cc0f9..9a80534e58 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -133,12 +133,13 @@ void Output::setup(int flag) // will not write on last step of run unless multiple of every // set next_dump_any to smallest next_dump // if no dumps, set next_dump_any to last+1 so will not influence next - // dump custom/cfg may invoke computes so wrap with clear/add + // wrap dumps that invoke computes with clear/add if (ndump && update->restrict_output == 0) { for (int idump = 0; idump < ndump; idump++) { if (strcmp(dump[idump]->style,"custom") == 0 || - strcmp(dump[idump]->style,"cfg") == 0) + strcmp(dump[idump]->style,"cfg") == 0 || + strcmp(dump[idump]->style,"local") == 0) modify->clearstep_compute(); if ((ntimestep % dump_every[idump] == 0 && last_dump[idump] != ntimestep) || last_dump[idump] < 0) { @@ -148,7 +149,8 @@ void Output::setup(int flag) next_dump[idump] = (ntimestep/dump_every[idump])*dump_every[idump] + dump_every[idump]; if (strcmp(dump[idump]->style,"custom") == 0 || - strcmp(dump[idump]->style,"cfg") == 0) + strcmp(dump[idump]->style,"cfg") == 0 || + strcmp(dump[idump]->style,"local") == 0) modify->addstep_compute(next_dump[idump]); if (idump) next_dump_any = MYMIN(next_dump_any,next_dump[idump]); else next_dump_any = next_dump[0]; @@ -201,19 +203,21 @@ void Output::setup(int flag) void Output::write(int ntimestep) { // next_dump does not force output on last step of run - // dump custom/cfg may invoke computes so wrap with clear/add + // wrap dumps that invoke computes with clear/add if (next_dump_any == ntimestep) { for (int idump = 0; idump < ndump; idump++) { if (next_dump[idump] == ntimestep && last_dump[idump] != ntimestep) { if (strcmp(dump[idump]->style,"custom") == 0 || - strcmp(dump[idump]->style,"cfg") == 0) + strcmp(dump[idump]->style,"cfg") == 0 || + strcmp(dump[idump]->style,"local") == 0) modify->clearstep_compute(); dump[idump]->write(); last_dump[idump] = ntimestep; next_dump[idump] += dump_every[idump]; if (strcmp(dump[idump]->style,"custom") == 0 || - strcmp(dump[idump]->style,"cfg") == 0) + strcmp(dump[idump]->style,"cfg") == 0 || + strcmp(dump[idump]->style,"local") == 0) modify->addstep_compute(next_dump[idump]); } if (idump) next_dump_any = MYMIN(next_dump_any,next_dump[idump]); diff --git a/src/style.h b/src/style.h index c30212c9aa..ce4aa7da52 100644 --- a/src/style.h +++ b/src/style.h @@ -150,6 +150,7 @@ ComputeStyle(temp/sphere,ComputeTempSphere) #include "dump_cfg.h" #include "dump_custom.h" #include "dump_dcd.h" +#include "dump_local.h" #include "dump_xyz.h" #endif @@ -158,6 +159,7 @@ DumpStyle(atom,DumpAtom) DumpStyle(cfg,DumpCFG) DumpStyle(custom,DumpCustom) DumpStyle(dcd,DumpDCD) +DumpStyle(local,DumpLocal) DumpStyle(xyz,DumpXYZ) #endif diff --git a/src/update.cpp b/src/update.cpp index afb43199d1..87cb472b94 100644 --- a/src/update.cpp +++ b/src/update.cpp @@ -247,6 +247,7 @@ void Update::reset_timestep(int narg, char **arg) modify->compute[i]->invoked_scalar = -1; modify->compute[i]->invoked_vector = -1; modify->compute[i]->invoked_peratom = -1; + modify->compute[i]->invoked_local = -1; } for (int i = 0; i < modify->ncompute; i++)