diff --git a/src/fix_ave_time.cpp b/src/fix_ave_time.cpp index 91d2047545..8c81114908 100644 --- a/src/fix_ave_time.cpp +++ b/src/fix_ave_time.cpp @@ -22,18 +22,21 @@ #include "modify.h" #include "compute.h" #include "group.h" +#include "memory.h" #include "error.h" using namespace LAMMPS_NS; enum{COMPUTE,FIX}; +enum{SCALAR,VECTOR,BOTH}; +enum{ONE,RUNNING,WINDOW}; /* ---------------------------------------------------------------------- */ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { - if (narg != 10) error->all("Illegal fix ave/time command"); + if (narg < 8) error->all("Illegal fix ave/time command"); MPI_Comm_rank(world,&me); @@ -49,16 +52,53 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) : id = new char[n]; strcpy(id,arg[7]); - int flag = atoi(arg[8]); + // option defaults - if (strcmp(arg[9],"NULL") == 0) fp = NULL; - else if (me == 0) { - fp = fopen(arg[9],"w"); - if (fp == NULL) { - char str[128]; - sprintf(str,"Cannot open fix ave/time file %s",arg[9]); - error->one(str); - } + sflag = 1; + vflag = 0; + fp = NULL; + ave = ONE; + + // optional args + + int iarg = 8; + while (iarg < narg) { + if (strcmp(arg[iarg],"type") == 0) { + if (iarg+2 > narg) error->all("Illegal fix ave/time command"); + if (strcmp(arg[iarg+1],"scalar") == 0) { + sflag = 1; + vflag = 0; + } else if (strcmp(arg[iarg+1],"vector") == 0) { + sflag = 0; + vflag = 1; + } else if (strcmp(arg[iarg+1],"both") == 0) sflag = vflag = 1; + else error->all("Illegal fix ave/time command"); + iarg += 2; + } else if (strcmp(arg[iarg],"file") == 0) { + if (iarg+2 > narg) error->all("Illegal fix ave/time command"); + if (me == 0) { + fp = fopen(arg[iarg+1],"w"); + if (fp == NULL) { + char str[128]; + sprintf(str,"Cannot open fix ave/time file %s",arg[iarg+1]); + error->one(str); + } + } + iarg += 2; + } else if (strcmp(arg[iarg],"ave") == 0) { + if (iarg+2 > narg) error->all("Illegal fix ave/time command"); + if (strcmp(arg[iarg+1],"one") == 0) ave = ONE; + else if (strcmp(arg[iarg+1],"running") == 0) ave = RUNNING; + else if (strcmp(arg[iarg+1],"window") == 0) ave = WINDOW; + else error->all("Illegal fix ave/time command"); + if (ave == WINDOW) { + if (iarg+3 > narg) error->all("Illegal fix ave/time command"); + nwindow = atoi(arg[iarg+2]); + if (nwindow <= 0) error->all("Illegal fix ave/time command"); + } + iarg += 2; + if (ave == WINDOW) iarg++; + } else error->all("Illegal fix ave/time command"); } // setup and error check @@ -76,11 +116,6 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) : if (ifix < 0) error->all("Fix ID for fix ave/time does not exist"); } - if (flag < 0 || flag > 2) error->all("Illegal fix ave/time command"); - sflag = vflag = 0; - if (flag == 0 || flag == 2) sflag = 1; - if (flag == 1 || flag == 2) vflag = 1; - if (which == COMPUTE) { if (sflag && modify->compute[icompute]->scalar_flag == 0) error->all("Fix ave/time compute does not calculate a scalar"); @@ -121,15 +156,25 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) : fprintf(fp,"TimeStep Scalar-value Vector-values\n"); } - vector = NULL; + // allocate memory for averaging + + vector = vector_total = NULL; if (vflag) { if (which == COMPUTE) size_vector = modify->compute[icompute]->size_vector; else size_vector = modify->fix[ifix]->size_vector; vector = new double[size_vector]; + vector_total = new double[size_vector]; + for (int i = 0; i < size_vector; i++) vector_total[i] = 0.0; } + scalar_list = NULL; + vector_list = NULL; + if (sflag && ave == WINDOW) scalar_list = new double[nwindow]; + if (vflag && ave == WINDOW) + vector_list = memory->create_2d_double_array(nwindow,size_vector, + "fix ave/time:vector_list"); + // enable this fix to produce a global scalar and/or vector - // initialize values to 0.0 since thermo may call them on first step if (sflag) scalar_flag = 1; if (vflag) vector_flag = 1; @@ -137,15 +182,23 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) : if (which == COMPUTE) extensive = modify->compute[icompute]->extensive; else extensive = modify->fix[ifix]->extensive; - scalar = 0.0; - for (int i = 0; i < size_vector; i++) vector[i] = 0.0; - - // nvalid = next step on which end_of_step does something + // initializations irepeat = 0; + iwindow = window_limit = 0; + norm = 0; + + // nvalid = next step on which end_of_step does something + // can be this timestep if multiple of nfreq and nrepeat = 1 + // else backup from next multiple of nfreq + nvalid = (update->ntimestep/nfreq)*nfreq + nfreq; - nvalid -= (nrepeat-1)*nevery; - if (nvalid <= update->ntimestep) + if (nvalid-nfreq == update->ntimestep && nrepeat == 1) + nvalid = update->ntimestep; + else + nvalid -= (nrepeat-1)*nevery; + + if (nvalid < update->ntimestep) error->all("Fix ave/time cannot be started on this timestep"); } @@ -157,6 +210,9 @@ FixAveTime::~FixAveTime() if (fp && me == 0) fclose(fp); delete [] compute; delete [] vector; + delete [] vector_total; + delete [] scalar_list; + memory->destroy_2d_double_array(vector_list); } /* ---------------------------------------------------------------------- */ @@ -204,6 +260,15 @@ void FixAveTime::init() } } +/* ---------------------------------------------------------------------- + only does something if nvalid = current timestep +------------------------------------------------------------------------- */ + +void FixAveTime::setup() +{ + end_of_step(); +} + /* ---------------------------------------------------------------------- */ void FixAveTime::end_of_step() @@ -243,45 +308,89 @@ void FixAveTime::end_of_step() vector[i] += fix->compute_vector(i); } + // done if irepeat < nrepeat + irepeat++; nvalid += nevery; + if (irepeat < nrepeat) return; - // output the results + // average the final result for the Nfreq timestep // reset irepeat and nvalid - if (irepeat == nrepeat) { - double repeat = nrepeat; - if (sflag) scalar /= repeat; - if (vflag) for (i = 0; i < size_vector; i++) vector[i] /= repeat; + double repeat = nrepeat; + if (sflag) scalar /= repeat; + if (vflag) for (i = 0; i < size_vector; i++) vector[i] /= repeat; - if (fp && me == 0) { - fprintf(fp,"%d",update->ntimestep); - if (sflag) fprintf(fp," %g",scalar); - if (vflag) - for (i = 0; i < size_vector; i++) fprintf(fp," %g",vector[i]); - fprintf(fp,"\n"); - fflush(fp); + irepeat = 0; + nvalid = update->ntimestep+nfreq - (nrepeat-1)*nevery; + + // if ave = ONE, only single Nfreq timestep value is needed + // if ave = RUNNING, combine with all previous Nfreq timestep values + // if ave = WINDOW, comine with nwindow most recent Nfreq timestep values + + if (ave == ONE) { + if (sflag) scalar_total = scalar; + if (vflag) for (i = 0; i < size_vector; i++) vector_total[i] = vector[i]; + norm = 1; + + } else if (ave == RUNNING) { + if (sflag) scalar_total += scalar; + if (vflag) for (i = 0; i < size_vector; i++) vector_total[i] += vector[i]; + norm++; + + } else if (ave == WINDOW) { + if (sflag) { + scalar_total += scalar; + if (window_limit) scalar_total -= scalar_list[iwindow]; + scalar_list[iwindow] = scalar; + } + if (vflag) { + for (i = 0; i < size_vector; i++) { + vector_total[i] += vector[i]; + if (window_limit) vector_total[i] -= vector_list[iwindow][i]; + vector_list[iwindow][i] = vector[i]; + } } - irepeat = 0; - nvalid = update->ntimestep+nfreq - (nrepeat-1)*nevery; + iwindow++; + if (iwindow == nwindow) { + iwindow = 0; + window_limit = 1; + } + if (window_limit) norm = nwindow; + else norm = iwindow; + } + + // output result to file + + if (fp && me == 0) { + fprintf(fp,"%d",update->ntimestep); + if (sflag) fprintf(fp," %g",scalar_total/norm); + if (vflag) + for (i = 0; i < size_vector; i++) fprintf(fp," %g",vector_total[i]/norm); + fprintf(fp,"\n"); + fflush(fp); } } /* ---------------------------------------------------------------------- return scalar value + could be ONE, RUNNING, or WINDOW value ------------------------------------------------------------------------- */ double FixAveTime::compute_scalar() { - return scalar; + if (norm) return scalar_total/norm; + return 0.0; } /* ---------------------------------------------------------------------- return Nth vector value + could be ONE, RUNNING, or WINDOW value ------------------------------------------------------------------------- */ double FixAveTime::compute_vector(int n) { - return vector[n]; + if (norm) return vector_total[n]/norm; + else return 0.0; } diff --git a/src/fix_ave_time.h b/src/fix_ave_time.h index 22a1dd72d9..22f3fecbc3 100644 --- a/src/fix_ave_time.h +++ b/src/fix_ave_time.h @@ -25,6 +25,7 @@ class FixAveTime : public Fix { ~FixAveTime(); int setmask(); void init(); + void setup(); void end_of_step(); double compute_scalar(); double compute_vector(int); @@ -35,11 +36,17 @@ class FixAveTime : public Fix { char *id; FILE *fp; - int sflag,vflag,nsum; + int sflag,vflag,ave,nwindow,nsum; double scalar,*vector; int ncompute; class Compute **compute; class Fix *fix; + + int norm,iwindow,window_limit; + double scalar_total; + double *scalar_list; + double *vector_total; + double **vector_list; }; } diff --git a/src/thermo.cpp b/src/thermo.cpp index fff54b9280..1778ab76c4 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -45,7 +45,6 @@ using namespace LAMMPS_NS; // vol, lx, ly, lz, xlo, xhi, ylo, yhi, zlo, zhi // pxx, pyy, pzz, pxy, pxz, pyz // drot, grot (rotational KE for dipole and granular particles) -// tave, pave, eave, peave (time-averaged quantities) // customize a new thermo style by adding a DEFINE to this list @@ -80,9 +79,6 @@ Thermo::Thermo(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) lostflag = ERROR; lostbefore = 0; flushflag = 0; - nwindow = 10; - ncurrent_t = ncurrent_p = ncurrent_e = ncurrent_pe = -1; - npartial_t = npartial_p = npartial_e = npartial_pe = 0; // set style and corresponding lineflag // custom style builds its own line of keywords @@ -156,14 +152,6 @@ Thermo::Thermo(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) format_f_def = (char *) "%14.4f"; format_int_user = NULL; format_float_user = NULL; - - // average quantities - - tsum = psum = esum = pesum = 0.0; - tpast = new double[nwindow]; - ppast = new double[nwindow]; - epast = new double[nwindow]; - pepast = new double[nwindow]; } /* ---------------------------------------------------------------------- */ @@ -186,13 +174,6 @@ Thermo::~Thermo() delete [] format_int_user; delete [] format_float_user; - - // average arrays - - delete [] tpast; - delete [] ppast; - delete [] epast; - delete [] pepast; } /* ---------------------------------------------------------------------- */ @@ -536,23 +517,6 @@ void Thermo::modify_params(int narg, char **arg) } iarg += 3; - } else if (strcmp(arg[iarg],"window") == 0) { - if (iarg+2 > narg) error->all("Illegal thermo_modify command"); - nwindow = atoi(arg[iarg+1]); - if (nwindow <= 0) error->all("Illegal thermo_modify command"); - ncurrent_t = ncurrent_p = ncurrent_e = ncurrent_pe = -1; - npartial_t = npartial_p = npartial_e = npartial_pe = 0; - tsum = psum = esum = pesum = 0.0; - delete [] tpast; - delete [] ppast; - delete [] epast; - delete [] pepast; - tpast = new double[nwindow]; - ppast = new double[nwindow]; - epast = new double[nwindow]; - pepast = new double[nwindow]; - iarg += 2; - } else error->all("Illegal thermo_modify command"); } } @@ -746,21 +710,6 @@ void Thermo::parse_fields(char *str) addfield("RotKEgrn",&Thermo::compute_grot,FLOAT); index_grot = add_compute(id_grot,0); - } else if (strcmp(word,"tave") == 0) { - addfield("T_ave",&Thermo::compute_tave,FLOAT); - index_temp = add_compute(id_temp,0); - } else if (strcmp(word,"pave") == 0) { - addfield("P_ave",&Thermo::compute_pave,FLOAT); - index_temp = add_compute(id_temp,0); - index_press = add_compute(id_press,0); - } else if (strcmp(word,"eave") == 0) { - addfield("E_ave",&Thermo::compute_eave,FLOAT); - peflag = 1; - index_temp = add_compute(id_temp,0); - } else if (strcmp(word,"peave") == 0) { - addfield("PE_ave",&Thermo::compute_peave,FLOAT); - peflag = 1; - // compute value = c_ID, fix value = f_ID, variable value = v_ID // if no trailing [], then arg is set to 0, else arg is between [] // copy = at most 8 chars of ID to pass to addfield @@ -1004,10 +953,6 @@ int Thermo::evaluate_keyword(char *word, double *answer) else if (strcmp(word,"drot") == 0) compute_drot(); else if (strcmp(word,"grot") == 0) compute_grot(); - else if (strcmp(word,"tave") == 0) compute_tave(); - else if (strcmp(word,"pave") == 0) compute_pave(); - else if (strcmp(word,"eave") == 0) compute_eave(); - else if (strcmp(word,"peave") == 0) compute_peave(); else return 1; thermoflag = 1; @@ -1459,99 +1404,3 @@ void Thermo::compute_grot() else dvalue = rotate_gran->compute_scalar(); if (normflag) dvalue /= natoms; } - -/* ---------------------------------------------------------------------- */ - -void Thermo::compute_tave() -{ - compute_temp(); - - if (firststep == 0) { - if (npartial_t == 0) { - ncurrent_t = 0; - npartial_t = 1; - tsum = tpast[0] = dvalue; - } - dvalue = tsum/npartial_t; - } else { - ncurrent_t++; - if (ncurrent_t == nwindow) ncurrent_t = 0; - if (npartial_t == nwindow) tsum -= tpast[ncurrent_t]; - else npartial_t++; - tpast[ncurrent_t] = dvalue; - tsum += tpast[ncurrent_t]; - dvalue = tsum/npartial_t; - } -} - -/* ---------------------------------------------------------------------- */ - -void Thermo::compute_pave() -{ - compute_press(); - - if (firststep == 0) { - if (npartial_p == 0) { - ncurrent_p = 0; - npartial_p = 1; - psum = ppast[0] = dvalue; - } - dvalue = psum/npartial_p; - } else { - ncurrent_p++; - if (ncurrent_p == nwindow) ncurrent_p = 0; - if (npartial_p == nwindow) psum -= ppast[ncurrent_p]; - else npartial_p++; - ppast[ncurrent_p] = dvalue; - psum += ppast[ncurrent_p]; - dvalue = psum/npartial_p; - } -} - -/* ---------------------------------------------------------------------- */ - -void Thermo::compute_eave() -{ - compute_etotal(); - - if (firststep == 0) { - if (npartial_e == 0) { - ncurrent_e = 0; - npartial_e = 1; - esum = epast[0] = dvalue; - } - dvalue = esum/npartial_e; - } else { - ncurrent_e++; - if (ncurrent_e == nwindow) ncurrent_e = 0; - if (npartial_e == nwindow) esum -= epast[ncurrent_e]; - else npartial_e++; - epast[ncurrent_e] = dvalue; - esum += epast[ncurrent_e]; - dvalue = esum/npartial_e; - } -} - -/* ---------------------------------------------------------------------- */ - -void Thermo::compute_peave() -{ - compute_pe(); - - if (firststep == 0) { - if (npartial_pe == 0) { - ncurrent_pe = 0; - npartial_pe = 1; - pesum = pepast[0] = dvalue; - } - dvalue = pesum/npartial_pe; - } else { - ncurrent_pe++; - if (ncurrent_pe == nwindow) ncurrent_pe = 0; - if (npartial_pe == nwindow) pesum -= pepast[ncurrent_pe]; - else npartial_pe++; - pepast[ncurrent_pe] = dvalue; - pesum += pepast[ncurrent_pe]; - dvalue = pesum/npartial_pe; - } -} diff --git a/src/thermo.h b/src/thermo.h index e1bf0f22d7..3f40731321 100644 --- a/src/thermo.h +++ b/src/thermo.h @@ -88,12 +88,6 @@ class Thermo : protected Pointers { int nvariable; // # of variables evaulated by thermo char **id_variable; // list of variable names - int nwindow; // time averaged values - int npartial_t,ncurrent_t,npartial_p,ncurrent_p; - int npartial_e,ncurrent_e,npartial_pe,ncurrent_pe; - double tsum,psum,esum,pesum; - double *tpast,*ppast,*epast,*pepast; - // private methods void allocate(); @@ -159,11 +153,6 @@ class Thermo : protected Pointers { void compute_drot(); void compute_grot(); - - void compute_tave(); - void compute_pave(); - void compute_eave(); - void compute_peave(); }; }