ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

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

Browse Source
This commit is contained in:
sjplimp
2007-10-18 21:47:59 +00:00
parent d7241b274d
commit ef7aaac85f
4 changed files with 156 additions and 202 deletions
Download Patch File Download Diff File Expand all files Collapse all files

171
src/fix_ave_time.cpp
View File

@ -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,17 +52,54 @@ 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");
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[9]);
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;
if (nvalid-nfreq == update->ntimestep && nrepeat == 1)
nvalid = update->ntimestep;
else
nvalid -= (nrepeat-1)*nevery;
if (nvalid <= update->ntimestep)
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;
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];
}
}
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;
}

9
src/fix_ave_time.h
View File

@ -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;
};
}

151
src/thermo.cpp
View File

@ -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;
}
}

11
src/thermo.h
View File

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