From 1afdd3c011d54a2a828eeb9dd8829bfe40f69a07 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Tue, 7 Dec 2021 09:16:19 -0700 Subject: [PATCH 01/13] new output vars for dumps --- src/dump.cpp | 19 +++++++++++++++ src/output.cpp | 66 ++++++++++++++++++++++++++++++++++---------------- src/output.h | 6 +++-- 3 files changed, 68 insertions(+), 23 deletions(-) diff --git a/src/dump.cpp b/src/dump.cpp index 4c02aa3070..639dc2fc07 100644 --- a/src/dump.cpp +++ b/src/dump.cpp @@ -937,9 +937,28 @@ void Dump::modify_params(int narg, char **arg) n = utils::inumeric(FLERR,arg[iarg+1],false,lmp); if (n <= 0) error->all(FLERR,"Illegal dump_modify command"); } + output->time_dump[idump] = 0; output->every_dump[idump] = n; iarg += 2; + } else if (strcmp(arg[iarg],"delta") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command"); + int idump; + for (idump = 0; idump < output->ndump; idump++) + if (strcmp(id,output->dump[idump]->id) == 0) break; + double delta; + if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) { + delete [] output->var_dump[idump]; + output->var_dump[idump] = utils::strdup(&arg[iarg+1][2]); + delta = 0.0; + } else { + delta = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (delta <= 0.0) error->all(FLERR,"Illegal dump_modify command"); + } + output->time_dump[idump] = 1; + output->delta_dump[idump] = delta; + iarg += 2; + } else if (strcmp(arg[iarg],"fileper") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command"); if (!multiproc) diff --git a/src/output.cpp b/src/output.cpp index 9e29caf118..0c283efa55 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -59,7 +59,9 @@ Output::Output(LAMMPS *lmp) : Pointers(lmp) ndump = 0; max_dump = 0; + time_dump = nullptr; every_dump = nullptr; + delta_dump = nullptr; next_dump = nullptr; last_dump = nullptr; var_dump = nullptr; @@ -92,7 +94,9 @@ Output::~Output() if (thermo) delete thermo; delete [] var_thermo; + memory->destroy(time_dump); memory->destroy(every_dump); + memory->destroy(delta_dump); memory->destroy(next_dump); memory->destroy(last_dump); for (int i = 0; i < ndump; i++) delete [] var_dump[i]; @@ -126,12 +130,13 @@ void Output::init() for (int i = 0; i < ndump; i++) dump[i]->init(); for (int i = 0; i < ndump; i++) - if (every_dump[i] == 0) { + if ((time_dump[i] == 0 && every_dump[i] == 0) || + (time_dump[i] == 1 && delta_dump[i] == 0.0)} { ivar_dump[i] = input->variable->find(var_dump[i]); if (ivar_dump[i] < 0) - error->all(FLERR,"Variable name for dump every does not exist"); + error->all(FLERR,"Variable name for dump every or delta does not exist"); if (!input->variable->equalstyle(ivar_dump[i])) - error->all(FLERR,"Variable for dump every is invalid style"); + error->all(FLERR,"Variable for dump every or delta is invalid style"); } if (restart_flag_single && restart_every_single == 0) { @@ -176,7 +181,7 @@ void Output::setup(int memflag) if (ndump && update->restrict_output == 0) { for (int idump = 0; idump < ndump; idump++) { - if (dump[idump]->clearstep || every_dump[idump] == 0) + if (dump[idump]->clearstep || var_dump[idump]) modify->clearstep_compute(); writeflag = 0; if (every_dump[idump] && ntimestep % every_dump[idump] == 0 && @@ -187,6 +192,7 @@ void Output::setup(int memflag) dump[idump]->write(); last_dump[idump] = ntimestep; } + if (every_dump[idump]) next_dump[idump] = (ntimestep/every_dump[idump])*every_dump[idump] + every_dump[idump]; @@ -197,7 +203,8 @@ void Output::setup(int memflag) error->all(FLERR,"Dump every variable returned a bad timestep"); next_dump[idump] = nextdump; } - if (dump[idump]->clearstep || every_dump[idump] == 0) { + + if (dump[idump]->clearstep || var_dump[idump]) { if (writeflag) modify->addstep_compute(next_dump[idump]); else modify->addstep_compute_all(next_dump[idump]); } @@ -287,21 +294,35 @@ void Output::write(bigint ntimestep) if (next_dump_any == ntimestep) { for (int idump = 0; idump < ndump; idump++) { if (next_dump[idump] == ntimestep) { - if (dump[idump]->clearstep || every_dump[idump] == 0) + if (dump[idump]->clearstep || var_dump[idump]) modify->clearstep_compute(); if (last_dump[idump] != ntimestep) { dump[idump]->write(); last_dump[idump] = ntimestep; } - if (every_dump[idump]) next_dump[idump] += every_dump[idump]; - else { - bigint nextdump = static_cast - (input->variable->compute_equal(ivar_dump[idump])); - if (nextdump <= ntimestep) - error->all(FLERR,"Dump every variable returned a bad timestep"); - next_dump[idump] = nextdump; + + if (time_dump[idump] == 0) { + if (every_dump[idump]) next_dump[idump] += every_dump[idump]; + else { + bigint nextdump = static_cast + (input->variable->compute_equal(ivar_dump[idump])); + if (nextdump <= ntimestep) + error->all(FLERR,"Dump every variable returned a bad timestep"); + next_dump[idump] = nextdump; + } + } else { + if (delta_dump[idump] > 0.0) { + bigint nextstep = + next_dump[idump] += every_dump[idump]; + } else { + double nexttime = input->variable->compute_equal(ivar_dump[idump]); + if (nexttime <= current_time) // NOTE: what is current time + error->all(FLERR,"Dump delta variable returned a bad time"); + next_dump[idump] = nextdump; + } } - if (dump[idump]->clearstep || every_dump[idump] == 0) + + if (dump[idump]->clearstep || var_dump[idump]) modify->addstep_compute(next_dump[idump]); } if (idump) next_dump_any = MIN(next_dump_any,next_dump[idump]); @@ -547,7 +568,9 @@ void Output::add_dump(int narg, char **arg) max_dump += DELTA; dump = (Dump **) memory->srealloc(dump,max_dump*sizeof(Dump *),"output:dump"); + memory->grow(time_dump,max_dump,"output:time_dump"); memory->grow(every_dump,max_dump,"output:every_dump"); + memory->grow(delta_dump,max_dump,"output:delta_dump"); memory->grow(next_dump,max_dump,"output:next_dump"); memory->grow(last_dump,max_dump,"output:last_dump"); var_dump = (char **) @@ -555,13 +578,6 @@ void Output::add_dump(int narg, char **arg) memory->grow(ivar_dump,max_dump,"output:ivar_dump"); } - // initialize per-dump data to suitable default values - - every_dump[ndump] = 0; - last_dump[ndump] = -1; - var_dump[ndump] = nullptr; - ivar_dump[ndump] = -1; - // create the Dump if (dump_map->find(arg[2]) != dump_map->end()) { @@ -569,10 +585,16 @@ void Output::add_dump(int narg, char **arg) dump[ndump] = dump_creator(lmp, narg, arg); } else error->all(FLERR,utils::check_packages_for_style("dump",arg[2],lmp)); + // initialize per-dump data to suitable default values + + time_dump[ndump] = 0; every_dump[ndump] = utils::inumeric(FLERR,arg[3],false,lmp); if (every_dump[ndump] <= 0) error->all(FLERR,"Illegal dump command"); + delta_dump[ndump] = 0.0; last_dump[ndump] = -1; var_dump[ndump] = nullptr; + ivar_dump[ndump] = -1; + ndump++; } @@ -624,7 +646,9 @@ void Output::delete_dump(char *id) for (int i = idump+1; i < ndump; i++) { dump[i-1] = dump[i]; + time_dump[i-1] = time_dump[i]; every_dump[i-1] = every_dump[i]; + delta_dump[i-1] = delta_dump[i]; next_dump[i-1] = next_dump[i]; last_dump[i-1] = last_dump[i]; var_dump[i-1] = var_dump[i]; diff --git a/src/output.h b/src/output.h index 9ae8b7fc3d..b47ec84dea 100644 --- a/src/output.h +++ b/src/output.h @@ -36,10 +36,12 @@ class Output : protected Pointers { int ndump; // # of Dumps defined int max_dump; // max size of Dump list bigint next_dump_any; // next timestep for any Dump - int *every_dump; // write freq for each Dump, 0 if var + int *time_dump; // 0/1 if write every N timesteps or Delta in sim time + int *every_dump; // write every this many timesteps, 0 if var + double *delta_dump; // write every this delta sim time, 0.0 if var bigint *next_dump; // next timestep to do each Dump bigint *last_dump; // last timestep each snapshot was output - char **var_dump; // variable name for dump frequency + char **var_dump; // variable name for dump frequency (steps or sim time) int *ivar_dump; // variable index for dump frequency Dump **dump; // list of defined Dumps From 26492b13d550f8b490e22e491c8b286f5ee5e4c4 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Tue, 7 Dec 2021 13:46:36 -0700 Subject: [PATCH 02/13] logic for dumps every steps and time delta --- src/dump.cpp | 42 ++++++------ src/output.cpp | 178 ++++++++++++++++++++++++++++++++++++------------- src/output.h | 14 ++-- 3 files changed, 160 insertions(+), 74 deletions(-) diff --git a/src/dump.cpp b/src/dump.cpp index 639dc2fc07..43fdb775b5 100644 --- a/src/dump.cpp +++ b/src/dump.cpp @@ -918,9 +918,22 @@ void Dump::modify_params(int narg, char **arg) else delay_flag = 0; iarg += 2; - } else if (strcmp(arg[iarg],"header") == 0) { + } else if (strcmp(arg[iarg],"delta") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command"); - header_flag = utils::logical(FLERR,arg[iarg+1],false,lmp); + int idump; + for (idump = 0; idump < output->ndump; idump++) + if (strcmp(id,output->dump[idump]->id) == 0) break; + double delta; + if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) { + delete [] output->var_dump[idump]; + output->var_dump[idump] = utils::strdup(&arg[iarg+1][2]); + delta = 0.0; + } else { + delta = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (delta <= 0.0) error->all(FLERR,"Illegal dump_modify command"); + } + output->mode_dump[idump] = 1; + output->delta_dump[idump] = delta; iarg += 2; } else if (strcmp(arg[iarg],"every") == 0) { @@ -937,28 +950,10 @@ void Dump::modify_params(int narg, char **arg) n = utils::inumeric(FLERR,arg[iarg+1],false,lmp); if (n <= 0) error->all(FLERR,"Illegal dump_modify command"); } - output->time_dump[idump] = 0; + output->mode_dump[idump] = 0; output->every_dump[idump] = n; iarg += 2; - } else if (strcmp(arg[iarg],"delta") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command"); - int idump; - for (idump = 0; idump < output->ndump; idump++) - if (strcmp(id,output->dump[idump]->id) == 0) break; - double delta; - if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) { - delete [] output->var_dump[idump]; - output->var_dump[idump] = utils::strdup(&arg[iarg+1][2]); - delta = 0.0; - } else { - delta = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (delta <= 0.0) error->all(FLERR,"Illegal dump_modify command"); - } - output->time_dump[idump] = 1; - output->delta_dump[idump] = delta; - iarg += 2; - } else if (strcmp(arg[iarg],"fileper") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command"); if (!multiproc) @@ -1026,6 +1021,11 @@ void Dump::modify_params(int narg, char **arg) iarg += n; } + } else if (strcmp(arg[iarg],"header") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command"); + header_flag = utils::logical(FLERR,arg[iarg+1],false,lmp); + iarg += 2; + } else if (strcmp(arg[iarg],"maxfiles") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command"); if (!multifile) diff --git a/src/output.cpp b/src/output.cpp index 0c283efa55..a47a4f8ec1 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -34,6 +34,7 @@ using namespace LAMMPS_NS; #define DELTA 1 +#define EPSDT 1.0e-6 /* ---------------------------------------------------------------------- initialize all output @@ -59,10 +60,11 @@ Output::Output(LAMMPS *lmp) : Pointers(lmp) ndump = 0; max_dump = 0; - time_dump = nullptr; + mode_dump = nullptr; every_dump = nullptr; delta_dump = nullptr; next_dump = nullptr; + next_time_dump = nullptr; last_dump = nullptr; var_dump = nullptr; ivar_dump = nullptr; @@ -94,10 +96,11 @@ Output::~Output() if (thermo) delete thermo; delete [] var_thermo; - memory->destroy(time_dump); + memory->destroy(mode_dump); memory->destroy(every_dump); memory->destroy(delta_dump); memory->destroy(next_dump); + memory->destroy(next_time_dump); memory->destroy(last_dump); for (int i = 0; i < ndump; i++) delete [] var_dump[i]; memory->sfree(var_dump); @@ -130,8 +133,8 @@ void Output::init() for (int i = 0; i < ndump; i++) dump[i]->init(); for (int i = 0; i < ndump; i++) - if ((time_dump[i] == 0 && every_dump[i] == 0) || - (time_dump[i] == 1 && delta_dump[i] == 0.0)} { + if ((mode_dump[i] == 0 && every_dump[i] == 0) || + (mode_dump[i] == 1 && delta_dump[i] == 0.0)) { ivar_dump[i] = input->variable->find(var_dump[i]); if (ivar_dump[i] < 0) error->all(FLERR,"Variable name for dump every or delta does not exist"); @@ -170,9 +173,9 @@ void Output::setup(int memflag) // current timestep is multiple of every and last dump not >= this step // this is first run after dump created and firstflag is set // note that variable freq will not write unless triggered by firstflag - // set next_dump to multiple of every or variable value + // set next_dump, and also next_time_dump for mode_dump = 1 // set next_dump_any to smallest next_dump - // wrap dumps that invoke computes and variable eval with clear/add + // wrap dumps that invoke computes or do variable eval with clear/add // if dump not written now, use addstep_compute_all() since don't know // what computes the dump write would invoke // if no dumps, set next_dump_any to last+1 so will not influence next @@ -183,9 +186,15 @@ void Output::setup(int memflag) for (int idump = 0; idump < ndump; idump++) { if (dump[idump]->clearstep || var_dump[idump]) modify->clearstep_compute(); + writeflag = 0; - if (every_dump[idump] && ntimestep % every_dump[idump] == 0 && - last_dump[idump] != ntimestep) writeflag = 1; + if (mode_dump[idump] == 0) { + if (every_dump[idump] && (ntimestep % every_dump[idump] == 0) && + last_dump[idump] != ntimestep) writeflag = 1; + } else { + if (delta_dump[idump] >= 0.0 && last_dump[idump] != ntimestep) + writeflag = 1; + } if (last_dump[idump] < 0 && dump[idump]->first_flag == 1) writeflag = 1; if (writeflag) { @@ -193,21 +202,15 @@ void Output::setup(int memflag) last_dump[idump] = ntimestep; } - if (every_dump[idump]) - next_dump[idump] = - (ntimestep/every_dump[idump])*every_dump[idump] + every_dump[idump]; - else { - bigint nextdump = static_cast - (input->variable->compute_equal(ivar_dump[idump])); - if (nextdump <= ntimestep) - error->all(FLERR,"Dump every variable returned a bad timestep"); - next_dump[idump] = nextdump; - } + // set next_dump and next_time_dump, 0 arg for setup() + + calculate_next_dump(0,idump,ntimestep); if (dump[idump]->clearstep || var_dump[idump]) { if (writeflag) modify->addstep_compute(next_dump[idump]); else modify->addstep_compute_all(next_dump[idump]); } + if (idump) next_dump_any = MIN(next_dump_any,next_dump[idump]); else next_dump_any = next_dump[0]; } @@ -288,43 +291,40 @@ void Output::setup(int memflag) void Output::write(bigint ntimestep) { - // next_dump does not force output on last step of run - // wrap dumps that invoke computes or eval of variable with clear/add + // perform dump if its next_dump = current ntimestep + // but not if it was already written on this step + // set next_dump and also next_time_dump for mode_dump = 1 + // set next_dump_any to smallest next_dump + // wrap dumps that invoke computes or do variable eval with clear/add + // if dump not written now, use addstep_compute_all() since don't know + // what computes the dump write would invoke + + int writeflag; if (next_dump_any == ntimestep) { for (int idump = 0; idump < ndump; idump++) { if (next_dump[idump] == ntimestep) { if (dump[idump]->clearstep || var_dump[idump]) modify->clearstep_compute(); - if (last_dump[idump] != ntimestep) { + + writeflag = 0; + if (last_dump[idump] != ntimestep) writeflag = 1; + + if (writeflag) { dump[idump]->write(); last_dump[idump] = ntimestep; } - if (time_dump[idump] == 0) { - if (every_dump[idump]) next_dump[idump] += every_dump[idump]; - else { - bigint nextdump = static_cast - (input->variable->compute_equal(ivar_dump[idump])); - if (nextdump <= ntimestep) - error->all(FLERR,"Dump every variable returned a bad timestep"); - next_dump[idump] = nextdump; - } - } else { - if (delta_dump[idump] > 0.0) { - bigint nextstep = - next_dump[idump] += every_dump[idump]; - } else { - double nexttime = input->variable->compute_equal(ivar_dump[idump]); - if (nexttime <= current_time) // NOTE: what is current time - error->all(FLERR,"Dump delta variable returned a bad time"); - next_dump[idump] = nextdump; - } - } + // set next_dump and next_time_dump, 1 arg for write() + + calculate_next_dump(1,idump,ntimestep); - if (dump[idump]->clearstep || var_dump[idump]) - modify->addstep_compute(next_dump[idump]); + if (dump[idump]->clearstep || var_dump[idump]) { + if (writeflag) modify->addstep_compute(next_dump[idump]); + else modify->addstep_compute_all(next_dump[idump]); + } } + if (idump) next_dump_any = MIN(next_dump_any,next_dump[idump]); else next_dump_any = next_dump[0]; } @@ -332,7 +332,7 @@ void Output::write(bigint ntimestep) // next_restart does not force output on last step of run // for toggle = 0, replace "*" with current timestep in restart filename - // eval of variable may invoke computes so wrap with clear/add + // next restart variable may invoke computes so wrap with clear/add if (next_restart == ntimestep) { if (next_restart_single == ntimestep) { @@ -417,6 +417,88 @@ void Output::write_dump(bigint ntimestep) } } +/* ---------------------------------------------------------------------- + calculate when next dump occurs for Dump instance idump + operates in one of two modes, based on mode_dump flag + for timestep mode, set next_dump + for simulation time mode, set next_time_dump and next_dump +------------------------------------------------------------------------- */ + + void Output::calculate_next_dump(int which, int idump, bigint ntimestep) +{ + // dump mode is by timestep + // just set next_dump + + if (mode_dump[idump] == 0) { + + // for setup, make next_dump a multiple of every_dump + + if (every_dump[idump]) { + if (which == 1) next_dump[idump] += every_dump[idump]; + else + next_dump[idump] = + (ntimestep/every_dump[idump])*every_dump[idump] + every_dump[idump]; + } else { + next_dump[idump] = static_cast + (input->variable->compute_equal(ivar_dump[idump])); + if (next_dump[idump] <= ntimestep) + error->all(FLERR,"Dump every variable returned a bad timestep"); + } + + // dump mode is by simulation time + // set next_time_dump and next_dump + + } else { + bigint nextdump; + double nexttime; + double tcurrent = update->atime + + (ntimestep - update->atimestep) * update->dt; + + // for setup, make nexttime a multiple of delta_dump + + if (delta_dump[idump] > 0.0) { + if (which == 1) nexttime = next_time_dump[idump] + delta_dump[idump]; + else + nexttime = static_cast (tcurrent/delta_dump[idump]) * + delta_dump[idump] + delta_dump[idump]; + + nextdump = ntimestep + + static_cast ((nexttime - tcurrent + EPSDT*update->dt) / + update->dt); + + // if delta is too small to reach next timestep, use multiple of delta + + if (nextdump == ntimestep) { + double tnext = update->atime + + (ntimestep+1 - update->atimestep) * update->dt; + int multiple = static_cast + ((tnext - nexttime) / delta_dump[idump]); + nexttime = nexttime + (multiple+1)*delta_dump[idump]; + nextdump = ntimestep + + static_cast ((nexttime - tcurrent + EPSDT*update->dt) / + update->dt); + } + + } else { + nexttime = input->variable->compute_equal(ivar_dump[idump]); + if (nexttime <= tcurrent) + error->all(FLERR,"Dump delta variable returned a bad time"); + nextdump = ntimestep + + static_cast ((nexttime - tcurrent + EPSDT*update->dt) / + update->dt); + if (nextdump <= ntimestep) + error->all(FLERR,"Dump delta variable too small for next timestep"); + } + + next_time_dump[idump] = nexttime; + next_dump[idump] = nextdump; + + //printf("END time %20.16g step %ld ratio %g\n", + // next_time_dump[idump],next_dump[idump], + // next_time_dump[idump]/update->dt/(next_dump[idump]+1)); + } +} + /* ---------------------------------------------------------------------- force restart file(s) to be written called from PRD and TAD @@ -568,10 +650,11 @@ void Output::add_dump(int narg, char **arg) max_dump += DELTA; dump = (Dump **) memory->srealloc(dump,max_dump*sizeof(Dump *),"output:dump"); - memory->grow(time_dump,max_dump,"output:time_dump"); + memory->grow(mode_dump,max_dump,"output:mode_dump"); memory->grow(every_dump,max_dump,"output:every_dump"); memory->grow(delta_dump,max_dump,"output:delta_dump"); memory->grow(next_dump,max_dump,"output:next_dump"); + memory->grow(next_time_dump,max_dump,"output:next_time_dump"); memory->grow(last_dump,max_dump,"output:last_dump"); var_dump = (char **) memory->srealloc(var_dump,max_dump*sizeof(char *),"output:var_dump"); @@ -587,7 +670,7 @@ void Output::add_dump(int narg, char **arg) // initialize per-dump data to suitable default values - time_dump[ndump] = 0; + mode_dump[ndump] = 0; every_dump[ndump] = utils::inumeric(FLERR,arg[3],false,lmp); if (every_dump[ndump] <= 0) error->all(FLERR,"Illegal dump command"); delta_dump[ndump] = 0.0; @@ -646,10 +729,11 @@ void Output::delete_dump(char *id) for (int i = idump+1; i < ndump; i++) { dump[i-1] = dump[i]; - time_dump[i-1] = time_dump[i]; + mode_dump[i-1] = mode_dump[i]; every_dump[i-1] = every_dump[i]; delta_dump[i-1] = delta_dump[i]; next_dump[i-1] = next_dump[i]; + next_time_dump[i-1] = next_time_dump[i]; last_dump[i-1] = last_dump[i]; var_dump[i-1] = var_dump[i]; ivar_dump[i-1] = ivar_dump[i]; diff --git a/src/output.h b/src/output.h index b47ec84dea..ed5bfc162b 100644 --- a/src/output.h +++ b/src/output.h @@ -36,13 +36,14 @@ class Output : protected Pointers { int ndump; // # of Dumps defined int max_dump; // max size of Dump list bigint next_dump_any; // next timestep for any Dump - int *time_dump; // 0/1 if write every N timesteps or Delta in sim time - int *every_dump; // write every this many timesteps, 0 if var - double *delta_dump; // write every this delta sim time, 0.0 if var - bigint *next_dump; // next timestep to do each Dump + int *mode_dump; // 0/1 if write every N timesteps or Delta in sim time + int *every_dump; // dump every this many timesteps, 0 if variable + double *delta_dump; // dump every this delta sim time, 0.0 if variable + bigint *next_dump; // next timestep to perform dump + double *next_time_dump; // next simulation time to perform dump (mode = 1) bigint *last_dump; // last timestep each snapshot was output - char **var_dump; // variable name for dump frequency (steps or sim time) - int *ivar_dump; // variable index for dump frequency + char **var_dump; // variable name for next dump (steps or sim time) + int *ivar_dump; // variable index of var_dump name Dump **dump; // list of defined Dumps int restart_flag; // 1 if any restart files are written @@ -89,6 +90,7 @@ class Output : protected Pointers { private: template static Dump *dump_creator(LAMMPS *, int, char **); + void calculate_next_dump(int, int, bigint); }; } // namespace LAMMPS_NS From d4149e913915399c51a2b93b61962847790539b0 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Wed, 8 Dec 2021 16:44:51 -0700 Subject: [PATCH 03/13] bug fixes to make a series of test inputs run correctly --- doc/src/dump.rst | 78 ++++++++------- doc/src/dump_modify.rst | 130 +++++++++++++++++++++---- doc/src/fix_dt_reset.rst | 18 +++- src/dump.cpp | 36 +++---- src/fix_dt_reset.cpp | 3 + src/output.cpp | 204 ++++++++++++++++++++++++++------------- src/output.h | 7 +- src/update.cpp | 26 +++-- 8 files changed, 342 insertions(+), 160 deletions(-) diff --git a/doc/src/dump.rst b/doc/src/dump.rst index c2509e6654..8c5ce6e208 100644 --- a/doc/src/dump.rst +++ b/doc/src/dump.rst @@ -169,11 +169,12 @@ or multiple smaller files). .. note:: - Because periodic boundary conditions are enforced only on - timesteps when neighbor lists are rebuilt, the coordinates of an atom - written to a dump file may be slightly outside the simulation box. - Re-neighbor timesteps will not typically coincide with the timesteps - dump snapshots are written. See the :doc:`dump_modify pbc ` command if you with to force coordinates to be + Because periodic boundary conditions are enforced only on timesteps + when neighbor lists are rebuilt, the coordinates of an atom written + to a dump file may be slightly outside the simulation box. + Re-neighbor timesteps will not typically coincide with the + timesteps dump snapshots are written. See the :doc:`dump_modify + pbc ` command if you with to force coordinates to be strictly inside the simulation box. .. note:: @@ -189,20 +190,21 @@ or multiple smaller files). multiple processors, each of which owns a subset of the atoms. For the *atom*, *custom*, *cfg*, and *local* styles, sorting is off by -default. For the *dcd*, *xtc*, *xyz*, and *molfile* styles, sorting by -atom ID is on by default. See the :doc:`dump_modify ` doc -page for details. +default. For the *dcd*, *xtc*, *xyz*, and *molfile* styles, sorting +by atom ID is on by default. See the :doc:`dump_modify ` +doc page for details. -The *atom/gz*, *cfg/gz*, *custom/gz*, *local/gz*, and *xyz/gz* styles are identical -in command syntax to the corresponding styles without "gz", however, -they generate compressed files using the zlib library. Thus the filename -suffix ".gz" is mandatory. This is an alternative approach to writing -compressed files via a pipe, as done by the regular dump styles, which -may be required on clusters where the interface to the high-speed network -disallows using the fork() library call (which is needed for a pipe). -For the remainder of this doc page, you should thus consider the *atom* -and *atom/gz* styles (etc) to be inter-changeable, with the exception -of the required filename suffix. +The *atom/gz*, *cfg/gz*, *custom/gz*, *local/gz*, and *xyz/gz* styles +are identical in command syntax to the corresponding styles without +"gz", however, they generate compressed files using the zlib +library. Thus the filename suffix ".gz" is mandatory. This is an +alternative approach to writing compressed files via a pipe, as done +by the regular dump styles, which may be required on clusters where +the interface to the high-speed network disallows using the fork() +library call (which is needed for a pipe). For the remainder of this +doc page, you should thus consider the *atom* and *atom/gz* styles +(etc) to be inter-changeable, with the exception of the required +filename suffix. Similarly, the *atom/zstd*, *cfg/zstd*, *custom/zstd*, *local/zstd*, and *xyz/zstd* styles are identical to the gz styles, but use the Zstd @@ -275,10 +277,11 @@ This bounding box is convenient for many visualization programs. The meaning of the 6 character flags for "xx yy zz" is the same as above. Note that the first two numbers on each line are now xlo_bound instead -of xlo, etc, since they represent a bounding box. See the :doc:`Howto triclinic ` page for a geometric description -of triclinic boxes, as defined by LAMMPS, simple formulas for how the -6 bounding box extents (xlo_bound,xhi_bound,etc) are calculated from -the triclinic parameters, and how to transform those parameters to and +of xlo, etc, since they represent a bounding box. See the :doc:`Howto +triclinic ` page for a geometric description of +triclinic boxes, as defined by LAMMPS, simple formulas for how the 6 +bounding box extents (xlo_bound,xhi_bound,etc) are calculated from the +triclinic parameters, and how to transform those parameters to and from other commonly used triclinic representations. The "ITEM: ATOMS" line in each snapshot lists column descriptors for @@ -310,23 +313,24 @@ written to the dump file. This local data is typically calculated by each processor based on the atoms it owns, but there may be zero or more entities per atom, e.g. a list of bond distances. An explanation of the possible dump local attributes is given below. Note that by -using input from the :doc:`compute property/local ` command with dump local, -it is possible to generate information on bonds, angles, etc that can -be cut and pasted directly into a data file read by the -:doc:`read_data ` command. +using input from the :doc:`compute property/local +` command with dump local, it is possible to +generate information on bonds, angles, etc that can be cut and pasted +directly into a data file read by the :doc:`read_data ` +command. Style *cfg* has the same command syntax as style *custom* and writes -extended CFG format files, as used by the -`AtomEye `_ visualization -package. Since the extended CFG format uses a single snapshot of the -system per file, a wildcard "\*" must be included in the filename, as -discussed below. The list of atom attributes for style *cfg* must -begin with either "mass type xs ys zs" or "mass type xsu ysu zsu" -since these quantities are needed to write the CFG files in the -appropriate format (though the "mass" and "type" fields do not appear -explicitly in the file). Any remaining attributes will be stored as -"auxiliary properties" in the CFG files. Note that you will typically -want to use the :doc:`dump_modify element ` command with +extended CFG format files, as used by the `AtomEye +`_ visualization package. +Since the extended CFG format uses a single snapshot of the system per +file, a wildcard "\*" must be included in the filename, as discussed +below. The list of atom attributes for style *cfg* must begin with +either "mass type xs ys zs" or "mass type xsu ysu zsu" since these +quantities are needed to write the CFG files in the appropriate format +(though the "mass" and "type" fields do not appear explicitly in the +file). Any remaining attributes will be stored as "auxiliary +properties" in the CFG files. Note that you will typically want to +use the :doc:`dump_modify element ` command with CFG-formatted files, to associate element names with atom types, so that AtomEye can render atoms appropriately. When unwrapped coordinates *xsu*, *ysu*, and *zsu* are requested, the nominal AtomEye diff --git a/doc/src/dump_modify.rst b/doc/src/dump_modify.rst index da7ccffeb2..42bb4d94fd 100644 --- a/doc/src/dump_modify.rst +++ b/doc/src/dump_modify.rst @@ -17,7 +17,7 @@ Syntax * one or more keyword/value pairs may be appended * these keywords apply to various dump styles -* keyword = *append* or *at* or *buffer* or *delay* or *element* or *every* or *fileper* or *first* or *flush* or *format* or *header* or *image* or *label* or *maxfiles* or *nfile* or *pad* or *pbc* or *precision* or *region* or *refresh* or *scale* or *sfactor* or *sort* or *tfactor* or *thermo* or *thresh* or *time* or *units* or *unwrap* +* keyword = *append* or *at* or *buffer* or *delay* or *element* or *every* or *every/time* or *fileper* or *first* or *flush* or *format* or *header* or *image* or *label* or *maxfiles* or *nfile* or *pad* or *pbc* or *precision* or *region* or *refresh* or *scale* or *sfactor* or *sort* or *tfactor* or *thermo* or *thresh* or *time* or *units* or *unwrap* .. parsed-literal:: @@ -32,6 +32,9 @@ Syntax *every* arg = N N = dump every this many timesteps N can be a variable (see below) + *every/time* arg = Delta + Delta = dump every this interval in simulation time (time units) + Delta can be a variable (see below) *fileper* arg = Np Np = write one file for every this many processors *first* arg = *yes* or *no* @@ -197,11 +200,18 @@ will be accepted. ---------- -The *every* keyword changes the dump frequency originally specified by -the :doc:`dump ` command to a new value. The every keyword can be -specified in one of two ways. It can be a numeric value in which case -it must be > 0. Or it can be an :doc:`equal-style variable `, -which should be specified as v_name, where name is the variable name. +The *every* keyword does two things. It specifies that the interval +between dump snapshots will be specified in timesteps, which is the +default if the *every* or *every/time* keywords are not used. See the +*every/time* keyword for how to specify the interval in simulation +time, i.e. in time units of the :doc:`units ` command. This +command also sets the interval value, which overrides the dump +frequency originally specified by the :doc:`dump ` command. + +The every keyword can be specified in one of two ways. It can be a +numeric value in which case it must be > 0. Or it can be an +:doc:`equal-style variable `, which should be specified as +v_name, where name is the variable name. In this case, the variable is evaluated at the beginning of a run to determine the next timestep at which a dump snapshot will be written @@ -210,11 +220,12 @@ determine the next timestep, etc. Thus the variable should return timestep values. See the stagger() and logfreq() and stride() math functions for :doc:`equal-style variables `, as examples of useful functions to use in this context. Other similar math functions -could easily be added as options for :doc:`equal-style variables `. Also see the next() function, which allows -use of a file-style variable which reads successive values from a -file, each time the variable is evaluated. Used with the *every* -keyword, if the file contains a list of ascending timesteps, you can -output snapshots whenever you wish. +could easily be added as options for :doc:`equal-style variables +`. Also see the next() function, which allows use of a +file-style variable which reads successive values from a file, each +time the variable is evaluated. Used with the *every* keyword, if the +file contains a list of ascending timesteps, you can output snapshots +whenever you wish. Note that when using the variable option with the *every* keyword, you need to use the *first* option if you want an initial snapshot written @@ -255,14 +266,95 @@ in file tmp.times: ---------- +The *every/time* keyword does two things. It specifies that the +interval between dump snapshots will be specified in simulation time, +i.e. in time units of the :doc:`units ` command. This can be +useful when the timestep size varies during a simulation run, e.g. by +use of the :doc:`fix dt/reset ` command. The default is +to specify the interval in timesteps; see the *every* keyword. This +command also sets the interval value. + +Note that since snapshots are output on simulation steps, each +snapshot will be written on the first timestep whose associated +simulation time is >= the exact snapshot time value. + +As with the *every* option, the *Delta* value can be specified in one +of two ways. It can be a numeric value in which case it must be > +0.0. Or it can be an :doc:`equal-style variable `, which +should be specified as v_name, where name is the variable name. + +In this case, the variable is evaluated at the beginning of a run to +determine the next simulation time at which a dump snapshot will be +written out. On that timestep the variable will be evaluated again to +determine the next simulation time, etc. Thus the variable should +return values in time units. Note the current timestep or simulation +time can be used in an :doc:`equal-style variables ` since +they are both thermodynamic keywords. Also see the next() function, +which allows use of a file-style variable which reads successive +values from a file, each time the variable is evaluated. Used with +the *every/time* keyword, if the file contains a list of ascending +simulation times, you can output snapshots whenever you wish. + +Note that when using the variable option with the *every/time* +keyword, you need to use the *first* option if you want an initial +snapshot written to the dump file. The *every/time* keyword cannot be +used with the dump *dcd* style. + +For example, the following commands will write snapshots at successive +simulation times which grow by a factor of 1.5 with each interval. +The dt value used in the variable is to avoid a zero result when the +initial simulation time is 0.0. + +.. code-block:: LAMMPS + + variable increase equal 1.5*(time+dt) + dump 1 all atom 100 tmp.dump + dump_modify 1 every/time v_increase first yes + +The following commands would write snapshots at the times listed in +file tmp.times: + +.. code-block:: LAMMPS + + variable f file tmp.times + variable s equal next(f) + dump 1 all atom 100 tmp.dump + dump_modify 1 every/time v_s + +.. note:: + + When using a file-style variable with the *every* keyword, the file + of timesteps must list a first times that is beyond the time + associated with the current timestep (e.g. it cannot be 0.0). And + it must list one or more times beyond the length of the run you + perform. This is because the dump command will generate an error + if the next time it reads from the file is not a value greater than + the current time. Thus if you wanted output at times 0,15,100 of a + run of length 100 in simulation time, the file should contain the + values 15,100,101 and you should also use the dump_modify first + command. Any final value > 100 could be used in place of 101. + +---------- + The *first* keyword determines whether a dump snapshot is written on the very first timestep after the dump command is invoked. This will -always occur if the current timestep is a multiple of N, the frequency -specified in the :doc:`dump ` command, including timestep 0. But -if this is not the case, a dump snapshot will only be written if the -setting of this keyword is *yes*\ . If it is *no*, which is the +always occur if the current timestep is a multiple of $N$, the +frequency specified in the :doc:`dump ` command or +:doc:`dump_modify every ` command, including timestep 0. +It will also always occur if the current simulation time is a multiple +of *Delta*, the time interval specified in the doc:`dump_modify every/time +` command. + +But if this is not the case, a dump snapshot will only be written if +the setting of this keyword is *yes*\ . If it is *no*, which is the default, then it will not be written. +Note that if the argument to the :doc:`dump_modify every +` doc:`dump_modify every/time ` commands is +a variable and not a numeric value, then specifying *first yes* is the +only way to write a dump snapshot on the first timestep after the dump +command is invoked. + ---------- The *flush* keyword determines whether a flush operation is invoked @@ -342,10 +434,10 @@ The *fileper* keyword is documented below with the *nfile* keyword. ---------- -The *header* keyword toggles whether the dump file will include a header. -Excluding a header will reduce the size of the dump file for fixes such as -:doc:`fix pair/tracker ` which do not require the information -typically written to the header. +The *header* keyword toggles whether the dump file will include a +header. Excluding a header will reduce the size of the dump file for +fixes such as :doc:`fix pair/tracker ` which do not +require the information typically written to the header. ---------- diff --git a/doc/src/fix_dt_reset.rst b/doc/src/fix_dt_reset.rst index c3aa431e18..368a3dcd70 100644 --- a/doc/src/fix_dt_reset.rst +++ b/doc/src/fix_dt_reset.rst @@ -78,13 +78,20 @@ outer loop (largest) timestep, which is the same timestep that the Note that the cumulative simulation time (in time units), which accounts for changes in the timestep size as a simulation proceeds, -can be accessed by the :doc:`thermo_style time ` keyword. +can be accessed by the :doc:`thermo_style time ` +keyword. + +Also note that the :doc:`dump_modify every/time ` option +allows dump files to be written at intervals specified by simulation +time, rather than by timesteps. Simulation time is in time units; +see the :doc:`units ` doc page for details. Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" -No information about this fix is written to :doc:`binary restart files `. None of the :doc:`fix_modify ` options -are relevant to this fix. +No information about this fix is written to :doc:`binary restart files +`. None of the :doc:`fix_modify ` options are +relevant to this fix. This fix computes a global scalar which can be accessed by various :doc:`output commands `. The scalar stores the last @@ -93,7 +100,8 @@ timestep on which the timestep was reset to a new value. The scalar value calculated by this fix is "intensive". No parameter of this fix can be used with the *start/stop* keywords of -the :doc:`run ` command. This fix is not invoked during :doc:`energy minimization `. +the :doc:`run ` command. This fix is not invoked during +:doc:`energy minimization `. Restrictions """""""""""" @@ -102,7 +110,7 @@ Restrictions Related commands """""""""""""""" -:doc:`timestep ` +:doc:`timestep `, :doc:`dump_modify every/time ` Default """"""" diff --git a/src/dump.cpp b/src/dump.cpp index 43fdb775b5..df39f3738d 100644 --- a/src/dump.cpp +++ b/src/dump.cpp @@ -918,24 +918,6 @@ void Dump::modify_params(int narg, char **arg) else delay_flag = 0; iarg += 2; - } else if (strcmp(arg[iarg],"delta") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command"); - int idump; - for (idump = 0; idump < output->ndump; idump++) - if (strcmp(id,output->dump[idump]->id) == 0) break; - double delta; - if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) { - delete [] output->var_dump[idump]; - output->var_dump[idump] = utils::strdup(&arg[iarg+1][2]); - delta = 0.0; - } else { - delta = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (delta <= 0.0) error->all(FLERR,"Illegal dump_modify command"); - } - output->mode_dump[idump] = 1; - output->delta_dump[idump] = delta; - iarg += 2; - } else if (strcmp(arg[iarg],"every") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command"); int idump; @@ -954,6 +936,24 @@ void Dump::modify_params(int narg, char **arg) output->every_dump[idump] = n; iarg += 2; + } else if (strcmp(arg[iarg],"every/time") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command"); + int idump; + for (idump = 0; idump < output->ndump; idump++) + if (strcmp(id,output->dump[idump]->id) == 0) break; + double delta; + if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) { + delete [] output->var_dump[idump]; + output->var_dump[idump] = utils::strdup(&arg[iarg+1][2]); + delta = 0.0; + } else { + delta = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (delta <= 0.0) error->all(FLERR,"Illegal dump_modify command"); + } + output->mode_dump[idump] = 1; + output->every_time_dump[idump] = delta; + iarg += 2; + } else if (strcmp(arg[iarg],"fileper") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command"); if (!multiproc) diff --git a/src/fix_dt_reset.cpp b/src/fix_dt_reset.cpp index c80c976504..fbbe7d7c3e 100644 --- a/src/fix_dt_reset.cpp +++ b/src/fix_dt_reset.cpp @@ -197,12 +197,15 @@ void FixDtReset::end_of_step() laststep = update->ntimestep; + // calls to other classes that need to know timestep size changed + update->update_time(); update->dt = dt; update->dt_default = 0; if (respaflag) update->integrate->reset_dt(); if (force->pair) force->pair->reset_dt(); for (int i = 0; i < modify->nfix; i++) modify->fix[i]->reset_dt(); + output->reset_dt(); } /* ---------------------------------------------------------------------- */ diff --git a/src/output.cpp b/src/output.cpp index a47a4f8ec1..9cc90ed4b5 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -29,6 +29,7 @@ #include "variable.h" #include "write_restart.h" +#include #include using namespace LAMMPS_NS; @@ -62,7 +63,7 @@ Output::Output(LAMMPS *lmp) : Pointers(lmp) max_dump = 0; mode_dump = nullptr; every_dump = nullptr; - delta_dump = nullptr; + every_time_dump = nullptr; next_dump = nullptr; next_time_dump = nullptr; last_dump = nullptr; @@ -98,7 +99,7 @@ Output::~Output() memory->destroy(mode_dump); memory->destroy(every_dump); - memory->destroy(delta_dump); + memory->destroy(every_time_dump); memory->destroy(next_dump); memory->destroy(next_time_dump); memory->destroy(last_dump); @@ -134,7 +135,7 @@ void Output::init() for (int i = 0; i < ndump; i++) dump[i]->init(); for (int i = 0; i < ndump; i++) if ((mode_dump[i] == 0 && every_dump[i] == 0) || - (mode_dump[i] == 1 && delta_dump[i] == 0.0)) { + (mode_dump[i] == 1 && every_time_dump[i] == 0.0)) { ivar_dump[i] = input->variable->find(var_dump[i]); if (ivar_dump[i] < 0) error->all(FLERR,"Variable name for dump every or delta does not exist"); @@ -169,42 +170,61 @@ void Output::setup(int memflag) { bigint ntimestep = update->ntimestep; - // perform dump at start of run only if: - // current timestep is multiple of every and last dump not >= this step - // this is first run after dump created and firstflag is set - // note that variable freq will not write unless triggered by firstflag - // set next_dump, and also next_time_dump for mode_dump = 1 - // set next_dump_any to smallest next_dump - // wrap dumps that invoke computes or do variable eval with clear/add - // if dump not written now, use addstep_compute_all() since don't know - // what computes the dump write would invoke - // if no dumps, set next_dump_any to last+1 so will not influence next - - int writeflag; + // consider all dumps + // decide whether to write snapshot and/or calculate next step for dump if (ndump && update->restrict_output == 0) { for (int idump = 0; idump < ndump; idump++) { + + // wrap dumps that invoke computes or do variable eval with clear/add + if (dump[idump]->clearstep || var_dump[idump]) modify->clearstep_compute(); - writeflag = 0; - if (mode_dump[idump] == 0) { - if (every_dump[idump] && (ntimestep % every_dump[idump] == 0) && - last_dump[idump] != ntimestep) writeflag = 1; - } else { - if (delta_dump[idump] >= 0.0 && last_dump[idump] != ntimestep) - writeflag = 1; - } + // write a snapshot at setup only if any of these 3 conditions hold + // (1) this is first run since dump was created and its first_flag = 0 + // (2) mode_dump = 0 and timestep is multiple of every_dump + // (3) mode_dump = 1 and time is multiple of every_time_dump (within EPSDT) + // (2) and (3) only apply for non-variable dump intervals + // finally, do not write if same snapshot written previously, + // i.e. on last timestep of previous run + + int writeflag = 0; + if (last_dump[idump] < 0 && dump[idump]->first_flag == 1) writeflag = 1; + if (mode_dump[idump] == 0) { + if (every_dump[idump] && (ntimestep % every_dump[idump] == 0)) + writeflag = 1; + } else { + if (every_time_dump[idump] > 0.0) { + double tcurrent = update->atime + + (ntimestep - update->atimestep) * update->dt; + double remainder = fmod(tcurrent,every_time_dump[idump]); + if ((remainder < EPSDT*update->dt) || + (every_time_dump[idump] - remainder < EPSDT*update->dt)) + writeflag = 1; + } + } + + if (last_dump[idump] == ntimestep) writeflag = 0; + + // perform dump + if (writeflag) { dump[idump]->write(); last_dump[idump] = ntimestep; } + // calculate timestep and/or time for next dump // set next_dump and next_time_dump, 0 arg for setup() + // only do this if dump written or dump has not been written yet - calculate_next_dump(0,idump,ntimestep); + if (writeflag || last_dump[idump] < 0) + calculate_next_dump(0,idump,ntimestep); + + // if dump not written now, use addstep_compute_all() + // since don't know what computes the dump will invoke if (dump[idump]->clearstep || var_dump[idump]) { if (writeflag) modify->addstep_compute(next_dump[idump]); @@ -214,6 +234,9 @@ void Output::setup(int memflag) if (idump) next_dump_any = MIN(next_dump_any,next_dump[idump]); else next_dump_any = next_dump[0]; } + + // if no dumps, set next_dump_any to last+1 so will not influence next + } else next_dump_any = update->laststep + 1; // do not write restart files at start of run @@ -296,33 +319,26 @@ void Output::write(bigint ntimestep) // set next_dump and also next_time_dump for mode_dump = 1 // set next_dump_any to smallest next_dump // wrap dumps that invoke computes or do variable eval with clear/add - // if dump not written now, use addstep_compute_all() since don't know - // what computes the dump write would invoke - + int writeflag; if (next_dump_any == ntimestep) { for (int idump = 0; idump < ndump; idump++) { if (next_dump[idump] == ntimestep) { + if (last_dump[idump] == ntimestep) continue; + if (dump[idump]->clearstep || var_dump[idump]) modify->clearstep_compute(); - writeflag = 0; - if (last_dump[idump] != ntimestep) writeflag = 1; + // perform dump + // reset next_dump and next_time_dump, 1 arg for write() - if (writeflag) { - dump[idump]->write(); - last_dump[idump] = ntimestep; - } - - // set next_dump and next_time_dump, 1 arg for write() - + dump[idump]->write(); + last_dump[idump] = ntimestep; calculate_next_dump(1,idump,ntimestep); - if (dump[idump]->clearstep || var_dump[idump]) { - if (writeflag) modify->addstep_compute(next_dump[idump]); - else modify->addstep_compute_all(next_dump[idump]); - } + if (dump[idump]->clearstep || var_dump[idump]) + modify->addstep_compute(next_dump[idump]); } if (idump) next_dump_any = MIN(next_dump_any,next_dump[idump]); @@ -422,6 +438,10 @@ void Output::write_dump(bigint ntimestep) operates in one of two modes, based on mode_dump flag for timestep mode, set next_dump for simulation time mode, set next_time_dump and next_dump + which flag depends on caller + 0 = from setup() at start of run + 1 = from write() during run each time a dump file is written + 2 = from reset_dt() called from fix dt/reset when it changes timestep size ------------------------------------------------------------------------- */ void Output::calculate_next_dump(int which, int idump, bigint ntimestep) @@ -431,13 +451,17 @@ void Output::write_dump(bigint ntimestep) if (mode_dump[idump] == 0) { - // for setup, make next_dump a multiple of every_dump - if (every_dump[idump]) { - if (which == 1) next_dump[idump] += every_dump[idump]; - else + + // which = 0: nextdump = next multiple of every_dump + // which = 1: increment nextdump by every_dump + + if (which == 0) next_dump[idump] = (ntimestep/every_dump[idump])*every_dump[idump] + every_dump[idump]; + else if (which == 1) + next_dump[idump] += every_dump[idump]; + } else { next_dump[idump] = static_cast (input->variable->compute_equal(ivar_dump[idump])); @@ -449,22 +473,29 @@ void Output::write_dump(bigint ntimestep) // set next_time_dump and next_dump } else { + bigint nextdump; double nexttime; double tcurrent = update->atime + (ntimestep - update->atimestep) * update->dt; - // for setup, make nexttime a multiple of delta_dump + if (every_time_dump[idump] > 0.0) { - if (delta_dump[idump] > 0.0) { - if (which == 1) nexttime = next_time_dump[idump] + delta_dump[idump]; - else - nexttime = static_cast (tcurrent/delta_dump[idump]) * - delta_dump[idump] + delta_dump[idump]; + // which = 0: nexttime = next multiple of every_time_dump + // which = 1: increment nexttime by every_time_dump + // which = 2: no change to previous nexttime (only timestep has changed) + + if (which == 0) + nexttime = static_cast (tcurrent/every_time_dump[idump]) * + every_time_dump[idump] + every_time_dump[idump]; + else if (which == 1) + nexttime = next_time_dump[idump] + every_time_dump[idump]; + else if (which == 2) + nexttime = next_time_dump[idump]; nextdump = ntimestep + - static_cast ((nexttime - tcurrent + EPSDT*update->dt) / - update->dt); + static_cast ((nexttime - tcurrent - EPSDT*update->dt) / + update->dt) + 1; // if delta is too small to reach next timestep, use multiple of delta @@ -472,30 +503,35 @@ void Output::write_dump(bigint ntimestep) double tnext = update->atime + (ntimestep+1 - update->atimestep) * update->dt; int multiple = static_cast - ((tnext - nexttime) / delta_dump[idump]); - nexttime = nexttime + (multiple+1)*delta_dump[idump]; + ((tnext - nexttime) / every_time_dump[idump]); + nexttime = nexttime + (multiple+1)*every_time_dump[idump]; nextdump = ntimestep + - static_cast ((nexttime - tcurrent + EPSDT*update->dt) / - update->dt); + static_cast ((nexttime - tcurrent - EPSDT*update->dt) / + update->dt) + 1; } } else { - nexttime = input->variable->compute_equal(ivar_dump[idump]); + + // do not re-evaulate variable for which = 2, leave nexttime as-is + // unless next_time_dump < 0.0, which means variable never yet evaluated + + if (which < 2 || next_time_dump[idump] < 0.0) + nexttime = input->variable->compute_equal(ivar_dump[idump]); + else + nexttime = next_time_dump[idump]; + if (nexttime <= tcurrent) - error->all(FLERR,"Dump delta variable returned a bad time"); + error->all(FLERR,"Dump every/time variable returned a bad time"); + nextdump = ntimestep + - static_cast ((nexttime - tcurrent + EPSDT*update->dt) / - update->dt); + static_cast ((nexttime - tcurrent - EPSDT*update->dt) / + update->dt) + 1; if (nextdump <= ntimestep) - error->all(FLERR,"Dump delta variable too small for next timestep"); + error->all(FLERR,"Dump every/time variable too small for next timestep"); } next_time_dump[idump] = nexttime; next_dump[idump] = nextdump; - - //printf("END time %20.16g step %ld ratio %g\n", - // next_time_dump[idump],next_dump[idump], - // next_time_dump[idump]/update->dt/(next_dump[idump]+1)); } } @@ -529,7 +565,7 @@ void Output::write_restart(bigint ntimestep) /* ---------------------------------------------------------------------- timestep is being changed, called by update->reset_timestep() - reset next timestep values for dumps, restart, thermo output + reset next output values for dumps, restart, thermo output reset to smallest value >= new timestep if next timestep set by variable evaluation, eval for ntimestep-1, so current ntimestep can be returned if needed @@ -626,6 +662,35 @@ void Output::reset_timestep(bigint ntimestep) next = MIN(next,next_thermo); } +/* ---------------------------------------------------------------------- + timestep size is being changed, called by fix dt/reset (at end of step) + reset next output values for dumps which have mode_dump=1 +------------------------------------------------------------------------- */ + +void Output::reset_dt() +{ + next_dump_any = MAXBIGINT; + + for (int idump = 0; idump < ndump; idump++) { + if (mode_dump[idump] == 1) { + + // reset next_dump for unchanged next_time_dump, 2 arg for reset_dt() + + calculate_next_dump(2,idump,update->ntimestep); + + // use compute_all() b/c don't know what computes will be needed + + if (dump[idump]->clearstep || var_dump[idump]) + modify->addstep_compute_all(next_dump[idump]); + } + + next_dump_any = MIN(next_dump_any,next_dump[idump]); + } + + next = MIN(next_dump_any,next_restart); + next = MIN(next,next_thermo); +} + /* ---------------------------------------------------------------------- add a Dump to list of Dumps ------------------------------------------------------------------------- */ @@ -652,7 +717,7 @@ void Output::add_dump(int narg, char **arg) memory->srealloc(dump,max_dump*sizeof(Dump *),"output:dump"); memory->grow(mode_dump,max_dump,"output:mode_dump"); memory->grow(every_dump,max_dump,"output:every_dump"); - memory->grow(delta_dump,max_dump,"output:delta_dump"); + memory->grow(every_time_dump,max_dump,"output:every_time_dump"); memory->grow(next_dump,max_dump,"output:next_dump"); memory->grow(next_time_dump,max_dump,"output:next_time_dump"); memory->grow(last_dump,max_dump,"output:last_dump"); @@ -673,7 +738,8 @@ void Output::add_dump(int narg, char **arg) mode_dump[ndump] = 0; every_dump[ndump] = utils::inumeric(FLERR,arg[3],false,lmp); if (every_dump[ndump] <= 0) error->all(FLERR,"Illegal dump command"); - delta_dump[ndump] = 0.0; + every_time_dump[ndump] = 0.0; + next_time_dump[ndump] = -1.0; last_dump[ndump] = -1; var_dump[ndump] = nullptr; ivar_dump[ndump] = -1; @@ -731,7 +797,7 @@ void Output::delete_dump(char *id) dump[i-1] = dump[i]; mode_dump[i-1] = mode_dump[i]; every_dump[i-1] = every_dump[i]; - delta_dump[i-1] = delta_dump[i]; + every_time_dump[i-1] = every_time_dump[i]; next_dump[i-1] = next_dump[i]; next_time_dump[i-1] = next_time_dump[i]; last_dump[i-1] = last_dump[i]; diff --git a/src/output.h b/src/output.h index ed5bfc162b..c8d3f734e0 100644 --- a/src/output.h +++ b/src/output.h @@ -37,8 +37,8 @@ class Output : protected Pointers { int max_dump; // max size of Dump list bigint next_dump_any; // next timestep for any Dump int *mode_dump; // 0/1 if write every N timesteps or Delta in sim time - int *every_dump; // dump every this many timesteps, 0 if variable - double *delta_dump; // dump every this delta sim time, 0.0 if variable + int *every_dump; // dump every N timesteps, 0 if variable + double *every_time_dump; // dump every Delta of sim time, 0.0 if variable bigint *next_dump; // next timestep to perform dump double *next_time_dump; // next simulation time to perform dump (mode = 1) bigint *last_dump; // last timestep each snapshot was output @@ -75,7 +75,8 @@ class Output : protected Pointers { void write(bigint); // output for current timestep void write_dump(bigint); // force output of dump snapshots void write_restart(bigint); // force output of a restart file - void reset_timestep(bigint); // reset next timestep for all output + void reset_timestep(bigint); // reset output which depeneds on timestep + void reset_dt(); // reset output which depends on dt void add_dump(int, char **); // add a Dump to Dump list void modify_dump(int, char **); // modify a Dump diff --git a/src/update.cpp b/src/update.cpp index 95dc47573e..10b4c573d7 100644 --- a/src/update.cpp +++ b/src/update.cpp @@ -458,7 +458,7 @@ Min *Update::minimize_creator(LAMMPS *lmp) } /* ---------------------------------------------------------------------- - reset timestep as called from input script + reset timestep called from input script ------------------------------------------------------------------------- */ void Update::reset_timestep(int narg, char **arg) @@ -470,24 +470,32 @@ void Update::reset_timestep(int narg, char **arg) /* ---------------------------------------------------------------------- reset timestep - called from rerun command and input script (indirectly) + called from input script (indirectly) or rerun command ------------------------------------------------------------------------- */ void Update::reset_timestep(bigint newstep) { + if (newstep < 0) error->all(FLERR,"Timestep must be >= 0"); + + bigint oldstep = ntimestep; ntimestep = newstep; - if (ntimestep < 0) error->all(FLERR,"Timestep must be >= 0"); - // set atimestep to new timestep - // so future update_time() calls will be correct + // if newstep >= oldstep, update simulation time accordingly + // if newstep < oldstep, zero simulation time - atimestep = ntimestep; + if (newstep >= oldstep) update_time(); - // trigger reset of timestep for output - // do not allow any timestep-dependent fixes to be already defined + if (newstep < oldstep) { + atime = 0.0; + atimestep = newstep; + } + + // changes to output that depend on timestep output->reset_timestep(ntimestep); + // do not allow timestep-dependent fixes to be defined + for (const auto &ifix : modify->get_fix_list()) if (ifix->time_depend) error->all(FLERR, "Cannot reset timestep with time-dependent fix {} defined",ifix->style); @@ -508,7 +516,7 @@ void Update::reset_timestep(bigint newstep) if (icompute->timeflag) icompute->clearstep(); } - // Neighbor Bin/Stencil/Pair classes store timestamps that need to be cleared + // neighbor Bin/Stencil/Pair classes store timestamps that need to be cleared neighbor->reset_timestep(ntimestep); } From e2969d09e1c892db45041dd3b619c80ee2554d1a Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Thu, 9 Dec 2021 11:53:47 -0700 Subject: [PATCH 04/13] bug fix for fix dt/reset freq of 1 --- src/output.cpp | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/src/output.cpp b/src/output.cpp index 9cc90ed4b5..fff3942780 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -515,9 +515,9 @@ void Output::write_dump(bigint ntimestep) // do not re-evaulate variable for which = 2, leave nexttime as-is // unless next_time_dump < 0.0, which means variable never yet evaluated - if (which < 2 || next_time_dump[idump] < 0.0) + if (which < 2 || next_time_dump[idump] < 0.0) { nexttime = input->variable->compute_equal(ivar_dump[idump]); - else + } else nexttime = next_time_dump[idump]; if (nexttime <= tcurrent) @@ -671,17 +671,20 @@ void Output::reset_dt() { next_dump_any = MAXBIGINT; + bigint ntimestep = update->ntimestep; + for (int idump = 0; idump < ndump; idump++) { if (mode_dump[idump] == 1) { // reset next_dump for unchanged next_time_dump, 2 arg for reset_dt() - - calculate_next_dump(2,idump,update->ntimestep); - + // do not invoke for a dump already scheduled for this step // use compute_all() b/c don't know what computes will be needed - if (dump[idump]->clearstep || var_dump[idump]) - modify->addstep_compute_all(next_dump[idump]); + if (next_dump[idump] != ntimestep) { + calculate_next_dump(2,idump,update->ntimestep); + if (dump[idump]->clearstep || var_dump[idump]) + modify->addstep_compute_all(next_dump[idump]); + } } next_dump_any = MIN(next_dump_any,next_dump[idump]); From 5ead32f886d58454a8cc9001c7119aabb702358c Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Fri, 10 Dec 2021 11:13:06 -0700 Subject: [PATCH 05/13] more debugging and features --- src/dump.cpp | 2 ++ src/fix_dt_reset.cpp | 1 + src/input.cpp | 16 +++++++++++++ src/output.cpp | 53 ++++++++++++++++++-------------------------- src/update.cpp | 1 + src/verlet.cpp | 2 ++ 6 files changed, 43 insertions(+), 32 deletions(-) diff --git a/src/dump.cpp b/src/dump.cpp index df39f3738d..21c1b16f75 100644 --- a/src/dump.cpp +++ b/src/dump.cpp @@ -330,6 +330,8 @@ void Dump::write() if (delay_flag && update->ntimestep < delaystep) return; + printf("DUMP %ld\n",update->ntimestep); + // if file per timestep, open new file if (multifile) openfile(); diff --git a/src/fix_dt_reset.cpp b/src/fix_dt_reset.cpp index fbbe7d7c3e..b7ab63f2d7 100644 --- a/src/fix_dt_reset.cpp +++ b/src/fix_dt_reset.cpp @@ -198,6 +198,7 @@ void FixDtReset::end_of_step() laststep = update->ntimestep; // calls to other classes that need to know timestep size changed + // similar logic is in Input::timestep() update->update_time(); update->dt = dt; diff --git a/src/input.cpp b/src/input.cpp index 08f62bdb42..35b04946c0 100644 --- a/src/input.cpp +++ b/src/input.cpp @@ -1836,8 +1836,24 @@ void Input::timer_command() void Input::timestep() { if (narg != 1) error->all(FLERR,"Illegal timestep command"); + + update->update_time(); update->dt = utils::numeric(FLERR,arg[0],false,lmp); update->dt_default = 0; + + if (update->first_update == 0) return; + + // calls to other classes that need to know timestep size changed + // similar logic is in FixDtReset::end_of_step() + // only if a run has already occurred + + int respaflag = 0; + if (utils::strmatch(update->integrate_style, "^respa")) respaflag = 1; + if (respaflag) update->integrate->reset_dt(); + + if (force->pair) force->pair->reset_dt(); + for (int i = 0; i < modify->nfix; i++) modify->fix[i]->reset_dt(); + output->reset_dt(); } /* ---------------------------------------------------------------------- */ diff --git a/src/output.cpp b/src/output.cpp index fff3942780..d34aa8d3f7 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -565,42 +565,22 @@ void Output::write_restart(bigint ntimestep) /* ---------------------------------------------------------------------- timestep is being changed, called by update->reset_timestep() - reset next output values for dumps, restart, thermo output + for dumps, require that no dump is "active" + meaning that a snapshot has already been output + reset next output values restart and thermo output reset to smallest value >= new timestep if next timestep set by variable evaluation, eval for ntimestep-1, so current ntimestep can be returned if needed no guarantee that variable can be evaluated for ntimestep-1 - if it depends on computes, but live with that rare case for now + e.g. if it depends on computes, but live with that rare case for now ------------------------------------------------------------------------- */ void Output::reset_timestep(bigint ntimestep) { next_dump_any = MAXBIGINT; - for (int idump = 0; idump < ndump; idump++) { - if (every_dump[idump]) { - next_dump[idump] = (ntimestep/every_dump[idump])*every_dump[idump]; - if (next_dump[idump] < ntimestep) next_dump[idump] += every_dump[idump]; - } else { - // ivar_dump may not be initialized - if (ivar_dump[idump] < 0) { - ivar_dump[idump] = input->variable->find(var_dump[idump]); - if (ivar_dump[idump] < 0) - error->all(FLERR,"Variable name for dump every does not exist"); - if (!input->variable->equalstyle(ivar_dump[idump])) - error->all(FLERR,"Variable for dump every is invalid style"); - } - modify->clearstep_compute(); - update->ntimestep--; - bigint nextdump = static_cast - (input->variable->compute_equal(ivar_dump[idump])); - if (nextdump < ntimestep) - error->all(FLERR,"Dump every variable returned a bad timestep"); - update->ntimestep++; - next_dump[idump] = nextdump; - modify->addstep_compute(next_dump[idump]); - } - next_dump_any = MIN(next_dump_any,next_dump[idump]); - } + for (int idump = 0; idump < ndump; idump++) + if (last_dump[idump] >= 0) + error->all("Cannot reset timestep with active dump - must undump first"); if (restart_flag_single) { if (restart_every_single) { @@ -669,21 +649,30 @@ void Output::reset_timestep(bigint ntimestep) void Output::reset_dt() { - next_dump_any = MAXBIGINT; - bigint ntimestep = update->ntimestep; + next_dump_any = MAXBIGINT; + for (int idump = 0; idump < ndump; idump++) { if (mode_dump[idump] == 1) { - // reset next_dump for unchanged next_time_dump, 2 arg for reset_dt() + // reset next_dump but do not change next_time_dump, 2 arg for reset_dt() // do not invoke for a dump already scheduled for this step // use compute_all() b/c don't know what computes will be needed if (next_dump[idump] != ntimestep) { calculate_next_dump(2,idump,update->ntimestep); - if (dump[idump]->clearstep || var_dump[idump]) - modify->addstep_compute_all(next_dump[idump]); + + // NOTE: think I should not do this here + // for time dumps, calc_next_dump should calc the next timestep + // as one less and not add it to computes + // then on that step, write() should not actually write the dump + // but trigger it on next step and addstep_compute_all for that step + // b/c when write() is called, the next-step timestep is set + // unless run every timestep is invoked in-between! + + //if (dump[idump]->clearstep || var_dump[idump]) + // modify->addstep_compute_all(next_dump[idump]); } } diff --git a/src/update.cpp b/src/update.cpp index 10b4c573d7..da412ee89e 100644 --- a/src/update.cpp +++ b/src/update.cpp @@ -491,6 +491,7 @@ void Update::reset_timestep(bigint newstep) } // changes to output that depend on timestep + // no active dumps allowed output->reset_timestep(ntimestep); diff --git a/src/verlet.cpp b/src/verlet.cpp index 342dc3d951..348639ef05 100644 --- a/src/verlet.cpp +++ b/src/verlet.cpp @@ -252,6 +252,8 @@ void Verlet::run(int n) ntimestep = ++update->ntimestep; ev_set(ntimestep); + printf("VERLET %ld: %d %d\n",ntimestep,eflag,vflag); + // initial time integration timer->stamp(); From 06c45fbe686f80d2b99442978031869618059a5f Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Mon, 20 Dec 2021 14:26:22 -0700 Subject: [PATCH 06/13] fix compiler errors --- src/input.cpp | 9 ++++++--- src/output.cpp | 3 ++- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/src/input.cpp b/src/input.cpp index 35b04946c0..30424ad5cb 100644 --- a/src/input.cpp +++ b/src/input.cpp @@ -27,9 +27,11 @@ #include "dihedral.h" #include "domain.h" #include "error.h" +#include "fix.h" #include "force.h" #include "group.h" #include "improper.h" +#include "integrate.h" #include "kspace.h" #include "memory.h" #include "min.h" @@ -1841,11 +1843,12 @@ void Input::timestep() update->dt = utils::numeric(FLERR,arg[0],false,lmp); update->dt_default = 0; - if (update->first_update == 0) return; - + // timestep command can be invoked between runs or by run every // calls to other classes that need to know timestep size changed // similar logic is in FixDtReset::end_of_step() - // only if a run has already occurred + // only need to do this if a run has already occurred + + if (update->first_update == 0) return; int respaflag = 0; if (utils::strmatch(update->integrate_style, "^respa")) respaflag = 1; diff --git a/src/output.cpp b/src/output.cpp index d34aa8d3f7..30248e449d 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -580,7 +580,8 @@ void Output::reset_timestep(bigint ntimestep) next_dump_any = MAXBIGINT; for (int idump = 0; idump < ndump; idump++) if (last_dump[idump] >= 0) - error->all("Cannot reset timestep with active dump - must undump first"); + error->all(FLERR, + "Cannot reset timestep with active dump - must undump first"); if (restart_flag_single) { if (restart_every_single) { From 4d31e300c61b8a2f72a10ee08136ad95ce0624ee Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Mon, 20 Dec 2021 16:39:17 -0700 Subject: [PATCH 07/13] change to checking timestep for time dumps at start of each step --- src/dump.cpp | 2 - src/integrate.cpp | 19 +- src/output.cpp | 824 +++++++++++++++++++++++----------------------- src/output.h | 6 +- src/verlet.cpp | 2 - 5 files changed, 440 insertions(+), 413 deletions(-) diff --git a/src/dump.cpp b/src/dump.cpp index 21c1b16f75..df39f3738d 100644 --- a/src/dump.cpp +++ b/src/dump.cpp @@ -330,8 +330,6 @@ void Dump::write() if (delay_flag && update->ntimestep < delaystep) return; - printf("DUMP %ld\n",update->ntimestep); - // if file per timestep, open new file if (multifile) openfile(); diff --git a/src/integrate.cpp b/src/integrate.cpp index 52d34fb943..7a1dbfaf34 100644 --- a/src/integrate.cpp +++ b/src/integrate.cpp @@ -20,6 +20,7 @@ #include "kspace.h" #include "modify.h" #include "pair.h" +#include "output.h" #include "update.h" using namespace LAMMPS_NS; @@ -115,7 +116,13 @@ void Integrate::ev_setup() /* ---------------------------------------------------------------------- set eflag,vflag for current iteration - based on computes that need energy/virial info on this timestep + based on + (1) computes that need energy/virial info on this timestep + (2) time dumps that need unknown per-atom info on this timestep + NOTE: could not check time dumps if timestep size is not varying + see NOTE in output.cpp + also inefficient to add all per-atom eng/virial computes + but don't know which ones the dump needs invoke matchstep() on all timestep-dependent computes to clear their arrays eflag: set any or no bits ENERGY_GLOBAL bit for global energy @@ -133,6 +140,10 @@ void Integrate::ev_set(bigint ntimestep) { int i,flag; + int tdflag = 0; + if (output->any_time_dumps) + tdflag = output->check_time_dumps(ntimestep); + flag = 0; int eflag_global = 0; for (i = 0; i < nelist_global; i++) @@ -143,7 +154,7 @@ void Integrate::ev_set(bigint ntimestep) int eflag_atom = 0; for (i = 0; i < nelist_atom; i++) if (elist_atom[i]->matchstep(ntimestep)) flag = 1; - if (flag) eflag_atom = ENERGY_ATOM; + if (flag || (tdflag && nelist_atom)) eflag_atom = ENERGY_ATOM; if (eflag_global) update->eflag_global = ntimestep; if (eflag_atom) update->eflag_atom = ntimestep; @@ -159,13 +170,13 @@ void Integrate::ev_set(bigint ntimestep) int vflag_atom = 0; for (i = 0; i < nvlist_atom; i++) if (vlist_atom[i]->matchstep(ntimestep)) flag = 1; - if (flag) vflag_atom = VIRIAL_ATOM; + if (flag || (tdflag && nvlist_atom)) vflag_atom = VIRIAL_ATOM; flag = 0; int cvflag_atom = 0; for (i = 0; i < ncvlist_atom; i++) if (cvlist_atom[i]->matchstep(ntimestep)) flag = 1; - if (flag) cvflag_atom = VIRIAL_CENTROID; + if (flag || (tdflag && ncvlist_atom)) cvflag_atom = VIRIAL_CENTROID; if (vflag_global) update->vflag_global = ntimestep; if (vflag_atom || cvflag_atom) update->vflag_atom = ntimestep; diff --git a/src/output.cpp b/src/output.cpp index 30248e449d..94688cd74c 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -133,7 +133,9 @@ void Output::init() } for (int i = 0; i < ndump; i++) dump[i]->init(); - for (int i = 0; i < ndump; i++) + any_time_dumps = 0; + for (int i = 0; i < ndump; i++) { + if (mode_dump[i]) any_time_dumps = 1; if ((mode_dump[i] == 0 && every_dump[i] == 0) || (mode_dump[i] == 1 && every_time_dump[i] == 0.0)) { ivar_dump[i] = input->variable->find(var_dump[i]); @@ -142,6 +144,7 @@ void Output::init() if (!input->variable->equalstyle(ivar_dump[i])) error->all(FLERR,"Variable for dump every or delta is invalid style"); } + } if (restart_flag_single && restart_every_single == 0) { ivar_restart_single = input->variable->find(var_restart_single); @@ -176,9 +179,11 @@ void Output::setup(int memflag) if (ndump && update->restrict_output == 0) { for (int idump = 0; idump < ndump; idump++) { - // wrap dumps that invoke computes or do variable eval with clear/add + // wrap step dumps that invoke computes or do variable eval with clear/add + // see NOTE in write() about also wrapping time dumps - if (dump[idump]->clearstep || var_dump[idump]) + if (mode_dump[idump] == 0 && + (dump[idump]->clearstep || var_dump[idump])) modify->clearstep_compute(); // write a snapshot at setup only if any of these 3 conditions hold @@ -226,7 +231,8 @@ void Output::setup(int memflag) // if dump not written now, use addstep_compute_all() // since don't know what computes the dump will invoke - if (dump[idump]->clearstep || var_dump[idump]) { + if (mode_dump[idump] == 0 && + (dump[idump]->clearstep || var_dump[idump])) { if (writeflag) modify->addstep_compute(next_dump[idump]); else modify->addstep_compute_all(next_dump[idump]); } @@ -239,413 +245,435 @@ void Output::setup(int memflag) } else next_dump_any = update->laststep + 1; - // do not write restart files at start of run - // set next_restart values to multiple of every or variable value - // wrap variable eval with clear/add - // if no restarts, set next_restart to last+1 so will not influence next + // do not write restart files at start of run + // set next_restart values to multiple of every or variable value + // wrap variable eval with clear/add + // if no restarts, set next_restart to last+1 so will not influence next - if (restart_flag && update->restrict_output == 0) { - if (restart_flag_single) { - if (restart_every_single) - next_restart_single = - (ntimestep/restart_every_single)*restart_every_single + - restart_every_single; - else { - bigint nextrestart = static_cast - (input->variable->compute_equal(ivar_restart_single)); - if (nextrestart <= ntimestep) - error->all(FLERR,"Restart variable returned a bad timestep"); - next_restart_single = nextrestart; - } - } else next_restart_single = update->laststep + 1; - if (restart_flag_double) { - if (restart_every_double) - next_restart_double = - (ntimestep/restart_every_double)*restart_every_double + - restart_every_double; - else { - bigint nextrestart = static_cast - (input->variable->compute_equal(ivar_restart_double)); - if (nextrestart <= ntimestep) - error->all(FLERR,"Restart variable returned a bad timestep"); - next_restart_double = nextrestart; - } - } else next_restart_double = update->laststep + 1; - next_restart = MIN(next_restart_single,next_restart_double); - } else next_restart = update->laststep + 1; + if (restart_flag && update->restrict_output == 0) { + if (restart_flag_single) { + if (restart_every_single) + next_restart_single = + (ntimestep/restart_every_single)*restart_every_single + + restart_every_single; + else { + bigint nextrestart = static_cast + (input->variable->compute_equal(ivar_restart_single)); + if (nextrestart <= ntimestep) + error->all(FLERR,"Restart variable returned a bad timestep"); + next_restart_single = nextrestart; + } + } else next_restart_single = update->laststep + 1; + if (restart_flag_double) { + if (restart_every_double) + next_restart_double = + (ntimestep/restart_every_double)*restart_every_double + + restart_every_double; + else { + bigint nextrestart = static_cast + (input->variable->compute_equal(ivar_restart_double)); + if (nextrestart <= ntimestep) + error->all(FLERR,"Restart variable returned a bad timestep"); + next_restart_double = nextrestart; + } + } else next_restart_double = update->laststep + 1; + next_restart = MIN(next_restart_single,next_restart_double); + } else next_restart = update->laststep + 1; - // print memory usage unless being called between multiple runs + // print memory usage unless being called between multiple runs - if (memflag) memory_usage(); + if (memflag) memory_usage(); - // set next_thermo to multiple of every or variable eval if var defined - // insure thermo output on last step of run - // thermo may invoke computes so wrap with clear/add + // set next_thermo to multiple of every or variable eval if var defined + // insure thermo output on last step of run + // thermo may invoke computes so wrap with clear/add - modify->clearstep_compute(); + modify->clearstep_compute(); - thermo->header(); - thermo->compute(0); - last_thermo = ntimestep; + thermo->header(); + thermo->compute(0); + last_thermo = ntimestep; - if (var_thermo) { - next_thermo = static_cast - (input->variable->compute_equal(ivar_thermo)); - if (next_thermo <= ntimestep) - error->all(FLERR,"Thermo every variable returned a bad timestep"); - } else if (thermo_every) { - next_thermo = (ntimestep/thermo_every)*thermo_every + thermo_every; - next_thermo = MIN(next_thermo,update->laststep); - } else next_thermo = update->laststep; + if (var_thermo) { + next_thermo = static_cast + (input->variable->compute_equal(ivar_thermo)); + if (next_thermo <= ntimestep) + error->all(FLERR,"Thermo every variable returned a bad timestep"); + } else if (thermo_every) { + next_thermo = (ntimestep/thermo_every)*thermo_every + thermo_every; + next_thermo = MIN(next_thermo,update->laststep); + } else next_thermo = update->laststep; - modify->addstep_compute(next_thermo); + modify->addstep_compute(next_thermo); - // next = next timestep any output will be done + // next = next timestep any output will be done - next = MIN(next_dump_any,next_restart); - next = MIN(next,next_thermo); -} + next = MIN(next_dump_any,next_restart); + next = MIN(next,next_thermo); + } -/* ---------------------------------------------------------------------- - perform all output for this timestep - only perform output if next matches current step and last output doesn't - do dump/restart before thermo so thermo CPU time will include them -------------------------------------------------------------------------- */ + /* ---------------------------------------------------------------------- + perform all output for this timestep + only perform output if next matches current step and last output doesn't + do dump/restart before thermo so thermo CPU time will include them + ------------------------------------------------------------------------- */ -void Output::write(bigint ntimestep) + void Output::write(bigint ntimestep) + { + // perform dump if its next_dump = current ntimestep + // but not if it was already written on this step + // set next_dump and also next_time_dump for mode_dump = 1 + // set next_dump_any to smallest next_dump + // wrap step dumps that invoke computes or do variable eval with clear/add + // NOTE: + // could wrap time dumps as well, if timestep size did not vary + // if wrap when timestep size varies frequently, + // then can do many unneeded addstep() --> inefficient + // hard to know if timestep varies, since run every could change it + // can't remove an uneeded addstep from a compute, b/c don't know + // what other command may have added it + + int writeflag; + + if (next_dump_any == ntimestep) { + for (int idump = 0; idump < ndump; idump++) { + if (next_dump[idump] == ntimestep) { + if (last_dump[idump] == ntimestep) continue; + + if (mode_dump[idump] == 0 && + (dump[idump]->clearstep || var_dump[idump])) + modify->clearstep_compute(); + + // perform dump + // reset next_dump and next_time_dump, 1 arg for write() + + dump[idump]->write(); + last_dump[idump] = ntimestep; + calculate_next_dump(1,idump,ntimestep); + + if (mode_dump[idump] == 0 && + (dump[idump]->clearstep || var_dump[idump])) + modify->addstep_compute(next_dump[idump]); + } + + if (idump) next_dump_any = MIN(next_dump_any,next_dump[idump]); + else next_dump_any = next_dump[0]; + } + } + + // next_restart does not force output on last step of run + // for toggle = 0, replace "*" with current timestep in restart filename + // next restart variable may invoke computes so wrap with clear/add + + if (next_restart == ntimestep) { + if (next_restart_single == ntimestep) { + + std::string file = restart1; + std::size_t found = file.find('*'); + if (found != std::string::npos) + file.replace(found,1,fmt::format("{}",update->ntimestep)); + + if (last_restart != ntimestep) restart->write(file); + + if (restart_every_single) next_restart_single += restart_every_single; + else { + modify->clearstep_compute(); + bigint nextrestart = static_cast + (input->variable->compute_equal(ivar_restart_single)); + if (nextrestart <= ntimestep) + error->all(FLERR,"Restart variable returned a bad timestep"); + next_restart_single = nextrestart; + modify->addstep_compute(next_restart_single); + } + } + if (next_restart_double == ntimestep) { + if (last_restart != ntimestep) { + if (restart_toggle == 0) { + restart->write(restart2a); + restart_toggle = 1; + } else { + restart->write(restart2b); + restart_toggle = 0; + } + } + if (restart_every_double) next_restart_double += restart_every_double; + else { + modify->clearstep_compute(); + bigint nextrestart = static_cast + (input->variable->compute_equal(ivar_restart_double)); + if (nextrestart <= ntimestep) + error->all(FLERR,"Restart variable returned a bad timestep"); + next_restart_double = nextrestart; + modify->addstep_compute(next_restart_double); + } + } + last_restart = ntimestep; + next_restart = MIN(next_restart_single,next_restart_double); + } + + // insure next_thermo forces output on last step of run + // thermo may invoke computes so wrap with clear/add + + if (next_thermo == ntimestep) { + modify->clearstep_compute(); + if (last_thermo != ntimestep) thermo->compute(1); + last_thermo = ntimestep; + if (var_thermo) { + next_thermo = static_cast + (input->variable->compute_equal(ivar_thermo)); + if (next_thermo <= ntimestep) + error->all(FLERR,"Thermo every variable returned a bad timestep"); + } else if (thermo_every) next_thermo += thermo_every; + else next_thermo = update->laststep; + next_thermo = MIN(next_thermo,update->laststep); + modify->addstep_compute(next_thermo); + } + + // next = next timestep any output will be done + + next = MIN(next_dump_any,next_restart); + next = MIN(next,next_thermo); + } + + /* ---------------------------------------------------------------------- + force a snapshot to be written for all dumps + called from PRD and TAD + ------------------------------------------------------------------------- */ + + void Output::write_dump(bigint ntimestep) + { + for (int idump = 0; idump < ndump; idump++) { + dump[idump]->write(); + last_dump[idump] = ntimestep; + } + } + + /* ---------------------------------------------------------------------- + calculate when next dump occurs for Dump instance idump + operates in one of two modes, based on mode_dump flag + for timestep mode, set next_dump + for simulation time mode, set next_time_dump and next_dump + which flag depends on caller + 0 = from setup() at start of run + 1 = from write() during run each time a dump file is written + 2 = from reset_dt() called from fix dt/reset when it changes timestep size + ------------------------------------------------------------------------- */ + + void Output::calculate_next_dump(int which, int idump, bigint ntimestep) + { + // dump mode is by timestep + // just set next_dump + + if (mode_dump[idump] == 0) { + + if (every_dump[idump]) { + + // which = 0: nextdump = next multiple of every_dump + // which = 1: increment nextdump by every_dump + + if (which == 0) + next_dump[idump] = + (ntimestep/every_dump[idump])*every_dump[idump] + every_dump[idump]; + else if (which == 1) + next_dump[idump] += every_dump[idump]; + + } else { + next_dump[idump] = static_cast + (input->variable->compute_equal(ivar_dump[idump])); + if (next_dump[idump] <= ntimestep) + error->all(FLERR,"Dump every variable returned a bad timestep"); + } + + // dump mode is by simulation time + // set next_time_dump and next_dump + + } else { + + bigint nextdump; + double nexttime; + double tcurrent = update->atime + + (ntimestep - update->atimestep) * update->dt; + + if (every_time_dump[idump] > 0.0) { + + // which = 0: nexttime = next multiple of every_time_dump + // which = 1: increment nexttime by every_time_dump + // which = 2: no change to previous nexttime (only timestep has changed) + + if (which == 0) + nexttime = static_cast (tcurrent/every_time_dump[idump]) * + every_time_dump[idump] + every_time_dump[idump]; + else if (which == 1) + nexttime = next_time_dump[idump] + every_time_dump[idump]; + else if (which == 2) + nexttime = next_time_dump[idump]; + + nextdump = ntimestep + + static_cast ((nexttime - tcurrent - EPSDT*update->dt) / + update->dt) + 1; + + // if delta is too small to reach next timestep, use multiple of delta + + if (nextdump == ntimestep) { + double tnext = update->atime + + (ntimestep+1 - update->atimestep) * update->dt; + int multiple = static_cast + ((tnext - nexttime) / every_time_dump[idump]); + nexttime = nexttime + (multiple+1)*every_time_dump[idump]; + nextdump = ntimestep + + static_cast ((nexttime - tcurrent - EPSDT*update->dt) / + update->dt) + 1; + } + + } else { + + // do not re-evaulate variable for which = 2, leave nexttime as-is + // unless next_time_dump < 0.0, which means variable never yet evaluated + + if (which < 2 || next_time_dump[idump] < 0.0) { + nexttime = input->variable->compute_equal(ivar_dump[idump]); + } else + nexttime = next_time_dump[idump]; + + if (nexttime <= tcurrent) + error->all(FLERR,"Dump every/time variable returned a bad time"); + + nextdump = ntimestep + + static_cast ((nexttime - tcurrent - EPSDT*update->dt) / + update->dt) + 1; + if (nextdump <= ntimestep) + error->all(FLERR,"Dump every/time variable too small for next timestep"); + } + + next_time_dump[idump] = nexttime; + next_dump[idump] = nextdump; + } + } + +/* ---------------------------------------------------------------------- */ + +int Output::check_time_dumps(bigint ntimestep) { - // perform dump if its next_dump = current ntimestep - // but not if it was already written on this step - // set next_dump and also next_time_dump for mode_dump = 1 - // set next_dump_any to smallest next_dump - // wrap dumps that invoke computes or do variable eval with clear/add + int nowflag = 0; + for (int i = 0; i < ndump; i++) + if (mode_dump[i] && next_dump[i] == ntimestep) nowflag = 1; - int writeflag; - - if (next_dump_any == ntimestep) { - for (int idump = 0; idump < ndump; idump++) { - if (next_dump[idump] == ntimestep) { - if (last_dump[idump] == ntimestep) continue; - - if (dump[idump]->clearstep || var_dump[idump]) - modify->clearstep_compute(); - - // perform dump - // reset next_dump and next_time_dump, 1 arg for write() - - dump[idump]->write(); - last_dump[idump] = ntimestep; - calculate_next_dump(1,idump,ntimestep); - - if (dump[idump]->clearstep || var_dump[idump]) - modify->addstep_compute(next_dump[idump]); - } - - if (idump) next_dump_any = MIN(next_dump_any,next_dump[idump]); - else next_dump_any = next_dump[0]; - } - } - - // next_restart does not force output on last step of run - // for toggle = 0, replace "*" with current timestep in restart filename - // next restart variable may invoke computes so wrap with clear/add - - if (next_restart == ntimestep) { - if (next_restart_single == ntimestep) { - - std::string file = restart1; - std::size_t found = file.find('*'); - if (found != std::string::npos) - file.replace(found,1,fmt::format("{}",update->ntimestep)); - - if (last_restart != ntimestep) restart->write(file); - - if (restart_every_single) next_restart_single += restart_every_single; - else { - modify->clearstep_compute(); - bigint nextrestart = static_cast - (input->variable->compute_equal(ivar_restart_single)); - if (nextrestart <= ntimestep) - error->all(FLERR,"Restart variable returned a bad timestep"); - next_restart_single = nextrestart; - modify->addstep_compute(next_restart_single); - } - } - if (next_restart_double == ntimestep) { - if (last_restart != ntimestep) { - if (restart_toggle == 0) { - restart->write(restart2a); - restart_toggle = 1; - } else { - restart->write(restart2b); - restart_toggle = 0; - } - } - if (restart_every_double) next_restart_double += restart_every_double; - else { - modify->clearstep_compute(); - bigint nextrestart = static_cast - (input->variable->compute_equal(ivar_restart_double)); - if (nextrestart <= ntimestep) - error->all(FLERR,"Restart variable returned a bad timestep"); - next_restart_double = nextrestart; - modify->addstep_compute(next_restart_double); - } - } - last_restart = ntimestep; - next_restart = MIN(next_restart_single,next_restart_double); - } - - // insure next_thermo forces output on last step of run - // thermo may invoke computes so wrap with clear/add - - if (next_thermo == ntimestep) { - modify->clearstep_compute(); - if (last_thermo != ntimestep) thermo->compute(1); - last_thermo = ntimestep; - if (var_thermo) { - next_thermo = static_cast - (input->variable->compute_equal(ivar_thermo)); - if (next_thermo <= ntimestep) - error->all(FLERR,"Thermo every variable returned a bad timestep"); - } else if (thermo_every) next_thermo += thermo_every; - else next_thermo = update->laststep; - next_thermo = MIN(next_thermo,update->laststep); - modify->addstep_compute(next_thermo); - } - - // next = next timestep any output will be done - - next = MIN(next_dump_any,next_restart); - next = MIN(next,next_thermo); + return nowflag; } -/* ---------------------------------------------------------------------- - force a snapshot to be written for all dumps - called from PRD and TAD -------------------------------------------------------------------------- */ + /* ---------------------------------------------------------------------- + force restart file(s) to be written + called from PRD and TAD + ------------------------------------------------------------------------- */ -void Output::write_dump(bigint ntimestep) -{ - for (int idump = 0; idump < ndump; idump++) { - dump[idump]->write(); - last_dump[idump] = ntimestep; - } -} + void Output::write_restart(bigint ntimestep) + { + if (restart_flag_single) { + std::string file = restart1; + std::size_t found = file.find('*'); + if (found != std::string::npos) + file.replace(found,1,fmt::format("{}",update->ntimestep)); + restart->write(file); + } + + if (restart_flag_double) { + if (restart_toggle == 0) { + restart->write(restart2a); + restart_toggle = 1; + } else { + restart->write(restart2b); + restart_toggle = 0; + } + } + + last_restart = ntimestep; + } + + /* ---------------------------------------------------------------------- + timestep is being changed, called by update->reset_timestep() + for dumps, require that no dump is "active" + meaning that a snapshot has already been output + reset next output values for restart and thermo + reset to smallest value >= new timestep + if next timestep set by variable evaluation, + eval for ntimestep-1, so current ntimestep can be returned if needed + no guarantee that variable can be evaluated for ntimestep-1 + e.g. if it depends on computes, but live with that rare case for now + ------------------------------------------------------------------------- */ + + void Output::reset_timestep(bigint ntimestep) + { + next_dump_any = MAXBIGINT; + for (int idump = 0; idump < ndump; idump++) + if (last_dump[idump] >= 0) + error->all(FLERR, + "Cannot reset timestep with active dump - must undump first"); + + if (restart_flag_single) { + if (restart_every_single) { + next_restart_single = + (ntimestep/restart_every_single)*restart_every_single; + if (next_restart_single < ntimestep) + next_restart_single += restart_every_single; + } else { + modify->clearstep_compute(); + update->ntimestep--; + bigint nextrestart = static_cast + (input->variable->compute_equal(ivar_restart_single)); + if (nextrestart < ntimestep) + error->all(FLERR,"Restart variable returned a bad timestep"); + update->ntimestep++; + next_restart_single = nextrestart; + modify->addstep_compute(next_restart_single); + } + } else next_restart_single = update->laststep + 1; + + if (restart_flag_double) { + if (restart_every_double) { + next_restart_double = + (ntimestep/restart_every_double)*restart_every_double; + if (next_restart_double < ntimestep) + next_restart_double += restart_every_double; + } else { + modify->clearstep_compute(); + update->ntimestep--; + bigint nextrestart = static_cast + (input->variable->compute_equal(ivar_restart_double)); + if (nextrestart < ntimestep) + error->all(FLERR,"Restart variable returned a bad timestep"); + update->ntimestep++; + next_restart_double = nextrestart; + modify->addstep_compute(next_restart_double); + } + } else next_restart_double = update->laststep + 1; + + next_restart = MIN(next_restart_single,next_restart_double); + + if (var_thermo) { + modify->clearstep_compute(); + update->ntimestep--; + next_thermo = static_cast + (input->variable->compute_equal(ivar_thermo)); + if (next_thermo < ntimestep) + error->all(FLERR,"Thermo_modify every variable returned a bad timestep"); + update->ntimestep++; + next_thermo = MIN(next_thermo,update->laststep); + modify->addstep_compute(next_thermo); + } else if (thermo_every) { + next_thermo = (ntimestep/thermo_every)*thermo_every; + if (next_thermo < ntimestep) next_thermo += thermo_every; + next_thermo = MIN(next_thermo,update->laststep); + } else next_thermo = update->laststep; + + next = MIN(next_dump_any,next_restart); + next = MIN(next,next_thermo); + } /* ---------------------------------------------------------------------- - calculate when next dump occurs for Dump instance idump - operates in one of two modes, based on mode_dump flag - for timestep mode, set next_dump - for simulation time mode, set next_time_dump and next_dump - which flag depends on caller - 0 = from setup() at start of run - 1 = from write() during run each time a dump file is written - 2 = from reset_dt() called from fix dt/reset when it changes timestep size -------------------------------------------------------------------------- */ - - void Output::calculate_next_dump(int which, int idump, bigint ntimestep) -{ - // dump mode is by timestep - // just set next_dump - - if (mode_dump[idump] == 0) { - - if (every_dump[idump]) { - - // which = 0: nextdump = next multiple of every_dump - // which = 1: increment nextdump by every_dump - - if (which == 0) - next_dump[idump] = - (ntimestep/every_dump[idump])*every_dump[idump] + every_dump[idump]; - else if (which == 1) - next_dump[idump] += every_dump[idump]; - - } else { - next_dump[idump] = static_cast - (input->variable->compute_equal(ivar_dump[idump])); - if (next_dump[idump] <= ntimestep) - error->all(FLERR,"Dump every variable returned a bad timestep"); - } - - // dump mode is by simulation time - // set next_time_dump and next_dump - - } else { - - bigint nextdump; - double nexttime; - double tcurrent = update->atime + - (ntimestep - update->atimestep) * update->dt; - - if (every_time_dump[idump] > 0.0) { - - // which = 0: nexttime = next multiple of every_time_dump - // which = 1: increment nexttime by every_time_dump - // which = 2: no change to previous nexttime (only timestep has changed) - - if (which == 0) - nexttime = static_cast (tcurrent/every_time_dump[idump]) * - every_time_dump[idump] + every_time_dump[idump]; - else if (which == 1) - nexttime = next_time_dump[idump] + every_time_dump[idump]; - else if (which == 2) - nexttime = next_time_dump[idump]; - - nextdump = ntimestep + - static_cast ((nexttime - tcurrent - EPSDT*update->dt) / - update->dt) + 1; - - // if delta is too small to reach next timestep, use multiple of delta - - if (nextdump == ntimestep) { - double tnext = update->atime + - (ntimestep+1 - update->atimestep) * update->dt; - int multiple = static_cast - ((tnext - nexttime) / every_time_dump[idump]); - nexttime = nexttime + (multiple+1)*every_time_dump[idump]; - nextdump = ntimestep + - static_cast ((nexttime - tcurrent - EPSDT*update->dt) / - update->dt) + 1; - } - - } else { - - // do not re-evaulate variable for which = 2, leave nexttime as-is - // unless next_time_dump < 0.0, which means variable never yet evaluated - - if (which < 2 || next_time_dump[idump] < 0.0) { - nexttime = input->variable->compute_equal(ivar_dump[idump]); - } else - nexttime = next_time_dump[idump]; - - if (nexttime <= tcurrent) - error->all(FLERR,"Dump every/time variable returned a bad time"); - - nextdump = ntimestep + - static_cast ((nexttime - tcurrent - EPSDT*update->dt) / - update->dt) + 1; - if (nextdump <= ntimestep) - error->all(FLERR,"Dump every/time variable too small for next timestep"); - } - - next_time_dump[idump] = nexttime; - next_dump[idump] = nextdump; - } -} - -/* ---------------------------------------------------------------------- - force restart file(s) to be written - called from PRD and TAD -------------------------------------------------------------------------- */ - -void Output::write_restart(bigint ntimestep) -{ - if (restart_flag_single) { - std::string file = restart1; - std::size_t found = file.find('*'); - if (found != std::string::npos) - file.replace(found,1,fmt::format("{}",update->ntimestep)); - restart->write(file); - } - - if (restart_flag_double) { - if (restart_toggle == 0) { - restart->write(restart2a); - restart_toggle = 1; - } else { - restart->write(restart2b); - restart_toggle = 0; - } - } - - last_restart = ntimestep; -} - -/* ---------------------------------------------------------------------- - timestep is being changed, called by update->reset_timestep() - for dumps, require that no dump is "active" - meaning that a snapshot has already been output - reset next output values restart and thermo output - reset to smallest value >= new timestep - if next timestep set by variable evaluation, - eval for ntimestep-1, so current ntimestep can be returned if needed - no guarantee that variable can be evaluated for ntimestep-1 - e.g. if it depends on computes, but live with that rare case for now -------------------------------------------------------------------------- */ - -void Output::reset_timestep(bigint ntimestep) -{ - next_dump_any = MAXBIGINT; - for (int idump = 0; idump < ndump; idump++) - if (last_dump[idump] >= 0) - error->all(FLERR, - "Cannot reset timestep with active dump - must undump first"); - - if (restart_flag_single) { - if (restart_every_single) { - next_restart_single = - (ntimestep/restart_every_single)*restart_every_single; - if (next_restart_single < ntimestep) - next_restart_single += restart_every_single; - } else { - modify->clearstep_compute(); - update->ntimestep--; - bigint nextrestart = static_cast - (input->variable->compute_equal(ivar_restart_single)); - if (nextrestart < ntimestep) - error->all(FLERR,"Restart variable returned a bad timestep"); - update->ntimestep++; - next_restart_single = nextrestart; - modify->addstep_compute(next_restart_single); - } - } else next_restart_single = update->laststep + 1; - - if (restart_flag_double) { - if (restart_every_double) { - next_restart_double = - (ntimestep/restart_every_double)*restart_every_double; - if (next_restart_double < ntimestep) - next_restart_double += restart_every_double; - } else { - modify->clearstep_compute(); - update->ntimestep--; - bigint nextrestart = static_cast - (input->variable->compute_equal(ivar_restart_double)); - if (nextrestart < ntimestep) - error->all(FLERR,"Restart variable returned a bad timestep"); - update->ntimestep++; - next_restart_double = nextrestart; - modify->addstep_compute(next_restart_double); - } - } else next_restart_double = update->laststep + 1; - - next_restart = MIN(next_restart_single,next_restart_double); - - if (var_thermo) { - modify->clearstep_compute(); - update->ntimestep--; - next_thermo = static_cast - (input->variable->compute_equal(ivar_thermo)); - if (next_thermo < ntimestep) - error->all(FLERR,"Thermo_modify every variable returned a bad timestep"); - update->ntimestep++; - next_thermo = MIN(next_thermo,update->laststep); - modify->addstep_compute(next_thermo); - } else if (thermo_every) { - next_thermo = (ntimestep/thermo_every)*thermo_every; - if (next_thermo < ntimestep) next_thermo += thermo_every; - next_thermo = MIN(next_thermo,update->laststep); - } else next_thermo = update->laststep; - - next = MIN(next_dump_any,next_restart); - next = MIN(next,next_thermo); -} - -/* ---------------------------------------------------------------------- - timestep size is being changed, called by fix dt/reset (at end of step) + timestep size is being changed reset next output values for dumps which have mode_dump=1 + called by fix dt/reset (at end of step) + or called by timestep command via run every (also at end of step) ------------------------------------------------------------------------- */ void Output::reset_dt() @@ -659,21 +687,10 @@ void Output::reset_dt() // reset next_dump but do not change next_time_dump, 2 arg for reset_dt() // do not invoke for a dump already scheduled for this step - // use compute_all() b/c don't know what computes will be needed + // since timestep change affects next step if (next_dump[idump] != ntimestep) { calculate_next_dump(2,idump,update->ntimestep); - - // NOTE: think I should not do this here - // for time dumps, calc_next_dump should calc the next timestep - // as one less and not add it to computes - // then on that step, write() should not actually write the dump - // but trigger it on next step and addstep_compute_all for that step - // b/c when write() is called, the next-step timestep is set - // unless run every timestep is invoked in-between! - - //if (dump[idump]->clearstep || var_dump[idump]) - // modify->addstep_compute_all(next_dump[idump]); } } @@ -684,6 +701,7 @@ void Output::reset_dt() next = MIN(next,next_thermo); } + /* ---------------------------------------------------------------------- add a Dump to list of Dumps ------------------------------------------------------------------------- */ diff --git a/src/output.h b/src/output.h index c8d3f734e0..c7c4faf252 100644 --- a/src/output.h +++ b/src/output.h @@ -36,6 +36,7 @@ class Output : protected Pointers { int ndump; // # of Dumps defined int max_dump; // max size of Dump list bigint next_dump_any; // next timestep for any Dump + int any_time_dumps; // 1 if any time dump defined int *mode_dump; // 0/1 if write every N timesteps or Delta in sim time int *every_dump; // dump every N timesteps, 0 if variable double *every_time_dump; // dump every Delta of sim time, 0.0 if variable @@ -75,13 +76,14 @@ class Output : protected Pointers { void write(bigint); // output for current timestep void write_dump(bigint); // force output of dump snapshots void write_restart(bigint); // force output of a restart file - void reset_timestep(bigint); // reset output which depeneds on timestep - void reset_dt(); // reset output which depends on dt + void reset_timestep(bigint); // reset output which depends on timestep + void reset_dt(); // reset output which depends on timestep size void add_dump(int, char **); // add a Dump to Dump list void modify_dump(int, char **); // modify a Dump void delete_dump(char *); // delete a Dump from Dump list int find_dump(const char *); // find a Dump ID + int check_time_dumps(bigint); // check if any time dump is output now void set_thermo(int, char **); // set thermo output freqquency void create_thermo(int, char **); // create a thermo style diff --git a/src/verlet.cpp b/src/verlet.cpp index 348639ef05..342dc3d951 100644 --- a/src/verlet.cpp +++ b/src/verlet.cpp @@ -252,8 +252,6 @@ void Verlet::run(int n) ntimestep = ++update->ntimestep; ev_set(ntimestep); - printf("VERLET %ld: %d %d\n",ntimestep,eflag,vflag); - // initial time integration timer->stamp(); From ded48cc031ffa255b2e17cc55f13ec04214dd4a3 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Tue, 21 Dec 2021 10:57:42 -0700 Subject: [PATCH 08/13] more optimizations and extend to other dump styles --- doc/src/dump_modify.rst | 49 +++++++++++++++++++++------------- src/COMPRESS/dump_xyz_gz.cpp | 8 +++++- src/COMPRESS/dump_xyz_zstd.cpp | 6 ++++- src/EXTRA-DUMP/dump_dcd.cpp | 8 ++++-- src/EXTRA-DUMP/dump_xtc.cpp | 25 ++++++++++------- src/dump.h | 4 +-- src/dump_xyz.cpp | 7 +++-- src/fix_dt_reset.cpp | 8 ------ src/integrate.cpp | 11 ++++---- src/output.cpp | 38 ++++++++++++++++++-------- src/output.h | 3 ++- 11 files changed, 106 insertions(+), 61 deletions(-) diff --git a/doc/src/dump_modify.rst b/doc/src/dump_modify.rst index 42bb4d94fd..0094fceb7b 100644 --- a/doc/src/dump_modify.rst +++ b/doc/src/dump_modify.rst @@ -200,15 +200,16 @@ will be accepted. ---------- -The *every* keyword does two things. It specifies that the interval -between dump snapshots will be specified in timesteps, which is the -default if the *every* or *every/time* keywords are not used. See the +The *every* keyword can be used with any dump style except the *dcd* +and *xtc* styles. It does two things. It specifies that the interval +between dump snapshots will be set in timesteps, which is the default +if the *every* or *every/time* keywords are not used. See the *every/time* keyword for how to specify the interval in simulation -time, i.e. in time units of the :doc:`units ` command. This -command also sets the interval value, which overrides the dump +time, i.e. in time units of the :doc:`units ` command. The +*every* keyword also sets the interval value, which overrides the dump frequency originally specified by the :doc:`dump ` command. -The every keyword can be specified in one of two ways. It can be a +The *every* keyword can be specified in one of two ways. It can be a numeric value in which case it must be > 0. Or it can be an :doc:`equal-style variable `, which should be specified as v_name, where name is the variable name. @@ -266,13 +267,21 @@ in file tmp.times: ---------- -The *every/time* keyword does two things. It specifies that the -interval between dump snapshots will be specified in simulation time, +The *every/time* keyword can be used with any dump style except the +*dcd* and *xtc* styles. It does two things. It specifies that the +interval between dump snapshots will be set in simulation time, i.e. in time units of the :doc:`units ` command. This can be useful when the timestep size varies during a simulation run, e.g. by use of the :doc:`fix dt/reset ` command. The default is -to specify the interval in timesteps; see the *every* keyword. This -command also sets the interval value. +to specify the interval in timesteps; see the *every* keyword. The +*every/time* command also sets the interval value. + +.. note:: + + If you wish dump styles *atom*, *custom*, *local*, or *xyz* to + include the simulation time as a field in the header portion of + each snapshot, you also need to use the dump_modify *time* keyword + with a setting of *yes*. See its documentation below. Note that since snapshots are output on simulation steps, each snapshot will be written on the first timestep whose associated @@ -323,8 +332,8 @@ file tmp.times: .. note:: - When using a file-style variable with the *every* keyword, the file - of timesteps must list a first times that is beyond the time + When using a file-style variable with the *every/time* keyword, the + file of timesteps must list a first time that is beyond the time associated with the current timestep (e.g. it cannot be 0.0). And it must list one or more times beyond the length of the run you perform. This is because the dump command will generate an error @@ -342,8 +351,8 @@ always occur if the current timestep is a multiple of $N$, the frequency specified in the :doc:`dump ` command or :doc:`dump_modify every ` command, including timestep 0. It will also always occur if the current simulation time is a multiple -of *Delta*, the time interval specified in the doc:`dump_modify every/time -` command. +of *Delta*, the time interval specified in the doc:`dump_modify +every/time ` command. But if this is not the case, a dump snapshot will only be written if the setting of this keyword is *yes*\ . If it is *no*, which is the @@ -731,16 +740,20 @@ threshold criterion is met. Otherwise it is not met. ---------- -The *time* keyword only applies to the dump *atom*, *custom*, and -*local* styles (and their COMPRESS package versions *atom/gz*, -*custom/gz* and *local/gz*\ ). If set to *yes*, each frame will will -contain two extra lines before the "ITEM: TIMESTEP" entry: +The *time* keyword only applies to the dump *atom*, *custom*, *local*, +and *xyz* styles (and their COMPRESS package versions *atom/gz*, +*custom/gz* and *local/gz*\ ). For the first 3 styles, if set to +*yes*, each frame will will contain two extra lines before the "ITEM: +TIMESTEP" entry: .. parsed-literal:: ITEM: TIME \ +For the *xyz* style, the simulation time is included on the same line +as the timestep value. + This will output the current elapsed simulation time in current time units equivalent to the :doc:`thermo keyword ` *time*\ . This is to simplify post-processing of trajectories using a variable time diff --git a/src/COMPRESS/dump_xyz_gz.cpp b/src/COMPRESS/dump_xyz_gz.cpp index 9d38c4673c..cec32090cb 100644 --- a/src/COMPRESS/dump_xyz_gz.cpp +++ b/src/COMPRESS/dump_xyz_gz.cpp @@ -88,11 +88,17 @@ void DumpXYZGZ::openfile() if (multifile) delete[] filecurrent; } +/* ---------------------------------------------------------------------- */ + void DumpXYZGZ::write_header(bigint ndump) { if (me == 0) { std::string header = fmt::format("{}\n", ndump); - header += fmt::format("Atoms. Timestep: {}\n", update->ntimestep); + if (time_flag) + header += fmt::format("Atoms. Timestep: {} Time: {:.6f}\n", + update->ntimestep, update->atime); + else + header += fmt::format("Atoms. Timestep: {}\n", update->ntimestep); writer.write(header.c_str(), header.length()); } } diff --git a/src/COMPRESS/dump_xyz_zstd.cpp b/src/COMPRESS/dump_xyz_zstd.cpp index bcbdc08a24..24d31fe671 100644 --- a/src/COMPRESS/dump_xyz_zstd.cpp +++ b/src/COMPRESS/dump_xyz_zstd.cpp @@ -100,7 +100,11 @@ void DumpXYZZstd::write_header(bigint ndump) { if (me == 0) { std::string header = fmt::format("{}\n", ndump); - header += fmt::format("Atoms. Timestep: {}\n", update->ntimestep); + if (time_flag) + header += fmt::format("Atoms. Timestep: {} Time: {:.6f}\n", + update->ntimestep, update->atime); + else + header += fmt::format("Atoms. Timestep: {}\n", update->ntimestep); writer.write(header.c_str(), header.length()); } } diff --git a/src/EXTRA-DUMP/dump_dcd.cpp b/src/EXTRA-DUMP/dump_dcd.cpp index ec3973448a..3aca5e8a98 100644 --- a/src/EXTRA-DUMP/dump_dcd.cpp +++ b/src/EXTRA-DUMP/dump_dcd.cpp @@ -100,13 +100,17 @@ void DumpDCD::init_style() if (sort_flag == 0 || sortcol != 0) error->all(FLERR,"Dump dcd requires sorting by atom ID"); - // check that dump frequency has not changed and is not a variable - // but only when not being called from the "write_dump" command. + // check that dump modify settings are compatible with dcd + // but only when not being called from the "write_dump" command if (strcmp(id,"WRITE_DUMP") != 0) { int idump; for (idump = 0; idump < output->ndump; idump++) if (strcmp(id,output->dump[idump]->id) == 0) break; + + if (output->mode_dump[idump] == 1) + error->all(FLERR,"Cannot use every/time setting for dump dcd"); + if (output->every_dump[idump] == 0) error->all(FLERR,"Cannot use variable every setting for dump dcd"); diff --git a/src/EXTRA-DUMP/dump_xtc.cpp b/src/EXTRA-DUMP/dump_xtc.cpp index 94846671cd..94271f31a6 100644 --- a/src/EXTRA-DUMP/dump_xtc.cpp +++ b/src/EXTRA-DUMP/dump_xtc.cpp @@ -121,17 +121,24 @@ void DumpXTC::init_style() if (flush_flag) error->all(FLERR,"Cannot set dump_modify flush for dump xtc"); - // check that dump frequency has not changed and is not a variable + // check that dump modify settings are compatible with xtc + // but only when not being called from the "write_dump" command - int idump; - for (idump = 0; idump < output->ndump; idump++) - if (strcmp(id,output->dump[idump]->id) == 0) break; - if (output->every_dump[idump] == 0) - error->all(FLERR,"Cannot use variable every setting for dump xtc"); + if (strcmp(id,"WRITE_DUMP") != 0) { + int idump; + for (idump = 0; idump < output->ndump; idump++) + if (strcmp(id,output->dump[idump]->id) == 0) break; + + if (output->mode_dump[idump] == 1) + error->all(FLERR,"Cannot use every/time setting for dump xtc"); - if (nevery_save == 0) nevery_save = output->every_dump[idump]; - else if (nevery_save != output->every_dump[idump]) - error->all(FLERR,"Cannot change dump_modify every for dump xtc"); + if (output->every_dump[idump] == 0) + error->all(FLERR,"Cannot use every variable setting for dump xtc"); + + if (nevery_save == 0) nevery_save = output->every_dump[idump]; + else if (nevery_save != output->every_dump[idump]) + error->all(FLERR,"Cannot change dump_modify every for dump xtc"); + } } /* ---------------------------------------------------------------------- */ diff --git a/src/dump.h b/src/dump.h index 730d4c9ca9..caeef827c7 100644 --- a/src/dump.h +++ b/src/dump.h @@ -26,7 +26,7 @@ class Dump : protected Pointers { int igroup, groupbit; // group that Dump is performed on int first_flag; // 0 if no initial dump, 1 if yes initial dump - int clearstep; // 1 if dump invokes computes, 0 if not + int clearstep; // 1 if dump can invoke computes, 0 if not int comm_forward; // size of forward communication (0 if none) int comm_reverse; // size of reverse communication (0 if none) @@ -75,7 +75,7 @@ class Dump : protected Pointers { int sortcol; // 0 to sort on ID, 1-N on columns int sortcolm1; // sortcol - 1 int sortorder; // ASCEND or DESCEND - int time_flag; // 1 if output accumulated time + int time_flag; // 1 if output simulation time int unit_flag; // 1 if dump should contain unit information int unit_count; // # of times the unit information was written int delay_flag; // 1 if delay output until delaystep diff --git a/src/dump_xyz.cpp b/src/dump_xyz.cpp index e009937959..db929f9d6e 100644 --- a/src/dump_xyz.cpp +++ b/src/dump_xyz.cpp @@ -131,7 +131,11 @@ void DumpXYZ::write_header(bigint n) { if (me == 0) { fprintf(fp,BIGINT_FORMAT "\n",n); - fprintf(fp,"Atoms. Timestep: " BIGINT_FORMAT "\n",update->ntimestep); + if (time_flag) + fprintf(fp,"Atoms. Timestep: " BIGINT_FORMAT " Time: %f\n", + update->ntimestep, update->atime); + else + fprintf(fp,"Atoms. Timestep: " BIGINT_FORMAT "\n",update->ntimestep); } } @@ -159,7 +163,6 @@ void DumpXYZ::pack(tagint *ids) } } - /* ---------------------------------------------------------------------- convert mybuf of doubles to one big formatted string in sbuf return -1 if strlen exceeds an int, since used as arg in MPI calls in Dump diff --git a/src/fix_dt_reset.cpp b/src/fix_dt_reset.cpp index b7ab63f2d7..adb0082fc8 100644 --- a/src/fix_dt_reset.cpp +++ b/src/fix_dt_reset.cpp @@ -121,14 +121,6 @@ void FixDtReset::init() respaflag = 0; if (utils::strmatch(update->integrate_style, "^respa")) respaflag = 1; - // check for DCD or XTC dumps - - for (int i = 0; i < output->ndump; i++) - if ((strcmp(output->dump[i]->style, "dcd") == 0 || - strcmp(output->dump[i]->style, "xtc") == 0) && - comm->me == 0) - error->warning(FLERR, "Dump dcd/xtc timestamp may be wrong with fix dt/reset"); - ftm2v = force->ftm2v; mvv2e = force->mvv2e; dt = update->dt; diff --git a/src/integrate.cpp b/src/integrate.cpp index 7a1dbfaf34..c111532ee2 100644 --- a/src/integrate.cpp +++ b/src/integrate.cpp @@ -118,11 +118,10 @@ void Integrate::ev_setup() set eflag,vflag for current iteration based on (1) computes that need energy/virial info on this timestep - (2) time dumps that need unknown per-atom info on this timestep - NOTE: could not check time dumps if timestep size is not varying + (2) time dumps that may need per-atom compute info on this timestep + NOTE: inefficient to add all per-atom eng/virial computes + but don't know which ones the dump needs see NOTE in output.cpp - also inefficient to add all per-atom eng/virial computes - but don't know which ones the dump needs invoke matchstep() on all timestep-dependent computes to clear their arrays eflag: set any or no bits ENERGY_GLOBAL bit for global energy @@ -141,8 +140,8 @@ void Integrate::ev_set(bigint ntimestep) int i,flag; int tdflag = 0; - if (output->any_time_dumps) - tdflag = output->check_time_dumps(ntimestep); + if (output->any_time_dumps && + output->next_time_dump_any == ntimestep) tdflag = 1; flag = 0; int eflag_global = 0; diff --git a/src/output.cpp b/src/output.cpp index 94688cd74c..d828b5a417 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -12,6 +12,10 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing author: Michal Kanski (Jagiellonian U) for simulation time dumps +------------------------------------------------------------------------- */ + #include "output.h" #include "style_dump.h" // IWYU pragma: keep @@ -177,6 +181,8 @@ void Output::setup(int memflag) // decide whether to write snapshot and/or calculate next step for dump if (ndump && update->restrict_output == 0) { + next_time_dump_any = MAXBIGINT; + for (int idump = 0; idump < ndump; idump++) { // wrap step dumps that invoke computes or do variable eval with clear/add @@ -237,6 +243,8 @@ void Output::setup(int memflag) else modify->addstep_compute_all(next_dump[idump]); } + if (mode_dump[idump] && (dump[idump]->clearstep || var_dump[idump])) + next_time_dump_any = MIN(next_time_dump_any,next_dump[idump]); if (idump) next_dump_any = MIN(next_dump_any,next_dump[idump]); else next_dump_any = next_dump[0]; } @@ -326,6 +334,9 @@ void Output::setup(int memflag) // set next_dump_any to smallest next_dump // wrap step dumps that invoke computes or do variable eval with clear/add // NOTE: + // not wrapping time dumps means that Integrate::ev_set() + // needs to trigger all per-atom eng/virial computes + // on a timestep where any time dump will be output // could wrap time dumps as well, if timestep size did not vary // if wrap when timestep size varies frequently, // then can do many unneeded addstep() --> inefficient @@ -336,7 +347,10 @@ void Output::setup(int memflag) int writeflag; if (next_dump_any == ntimestep) { + for (int idump = 0; idump < ndump; idump++) { + next_time_dump_any = MAXBIGINT; + if (next_dump[idump] == ntimestep) { if (last_dump[idump] == ntimestep) continue; @@ -356,6 +370,8 @@ void Output::setup(int memflag) modify->addstep_compute(next_dump[idump]); } + if (mode_dump[idump] && (dump[idump]->clearstep || var_dump[idump])) + next_time_dump_any = MIN(next_time_dump_any,next_dump[idump]); if (idump) next_dump_any = MIN(next_dump_any,next_dump[idump]); else next_dump_any = next_dump[0]; } @@ -680,23 +696,23 @@ void Output::reset_dt() { bigint ntimestep = update->ntimestep; - next_dump_any = MAXBIGINT; + next_time_dump_any = MAXBIGINT; for (int idump = 0; idump < ndump; idump++) { - if (mode_dump[idump] == 1) { + if (mode_dump[idump] == 0) continue; - // reset next_dump but do not change next_time_dump, 2 arg for reset_dt() - // do not invoke for a dump already scheduled for this step - // since timestep change affects next step + // reset next_dump but do not change next_time_dump, 2 arg for reset_dt() + // do not invoke for a dump already scheduled for this step + // since timestep change affects next step + + if (next_dump[idump] != ntimestep) + calculate_next_dump(2,idump,update->ntimestep); - if (next_dump[idump] != ntimestep) { - calculate_next_dump(2,idump,update->ntimestep); - } - } - - next_dump_any = MIN(next_dump_any,next_dump[idump]); + if (dump[idump]->clearstep || var_dump[idump]) + next_time_dump_any = MIN(next_time_dump_any,next_dump[idump]); } + next_dump_any = MIN(next_dump_any,next_time_dump_any); next = MIN(next_dump_any,next_restart); next = MIN(next,next_thermo); } diff --git a/src/output.h b/src/output.h index c7c4faf252..3f557bdef5 100644 --- a/src/output.h +++ b/src/output.h @@ -35,7 +35,8 @@ class Output : protected Pointers { int ndump; // # of Dumps defined int max_dump; // max size of Dump list - bigint next_dump_any; // next timestep for any Dump + bigint next_dump_any; // next timestep for any dump + bigint next_time_dump_any; // next timestep for any time dump with computes int any_time_dumps; // 1 if any time dump defined int *mode_dump; // 0/1 if write every N timesteps or Delta in sim time int *every_dump; // dump every N timesteps, 0 if variable From 576e7878399c02caeb426fca26be32ff27b32f3f Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Tue, 21 Dec 2021 12:18:26 -0700 Subject: [PATCH 09/13] make xyz dumps print out current simulation time --- src/COMPRESS/dump_xyz_gz.cpp | 11 +++++++---- src/COMPRESS/dump_xyz_zstd.cpp | 15 ++++++++++----- src/dump_xyz.cpp | 9 ++++++--- 3 files changed, 23 insertions(+), 12 deletions(-) diff --git a/src/COMPRESS/dump_xyz_gz.cpp b/src/COMPRESS/dump_xyz_gz.cpp index cec32090cb..a43a27ca7d 100644 --- a/src/COMPRESS/dump_xyz_gz.cpp +++ b/src/COMPRESS/dump_xyz_gz.cpp @@ -94,11 +94,14 @@ void DumpXYZGZ::write_header(bigint ndump) { if (me == 0) { std::string header = fmt::format("{}\n", ndump); - if (time_flag) - header += fmt::format("Atoms. Timestep: {} Time: {:.6f}\n", - update->ntimestep, update->atime); - else + if (time_flag) { + double tcurrent = update->atime + + (update->ntimestep-update->atimestep) + update->dt; + fprintf(fp,"Atoms. Timestep: " BIGINT_FORMAT " Time: %f\n", + update->ntimestep, tcurrent); + } else { header += fmt::format("Atoms. Timestep: {}\n", update->ntimestep); + } writer.write(header.c_str(), header.length()); } } diff --git a/src/COMPRESS/dump_xyz_zstd.cpp b/src/COMPRESS/dump_xyz_zstd.cpp index 24d31fe671..e955a0e3ba 100644 --- a/src/COMPRESS/dump_xyz_zstd.cpp +++ b/src/COMPRESS/dump_xyz_zstd.cpp @@ -96,15 +96,20 @@ void DumpXYZZstd::openfile() if (multifile) delete[] filecurrent; } +/* ---------------------------------------------------------------------- */ + void DumpXYZZstd::write_header(bigint ndump) { if (me == 0) { std::string header = fmt::format("{}\n", ndump); - if (time_flag) - header += fmt::format("Atoms. Timestep: {} Time: {:.6f}\n", - update->ntimestep, update->atime); - else - header += fmt::format("Atoms. Timestep: {}\n", update->ntimestep); + if (time_flag) { + double tcurrent = update->atime + + (update->ntimestep-update->atimestep) + update->dt; + fprintf(fp,"Atoms. Timestep: " BIGINT_FORMAT " Time: %f\n", + update->ntimestep, tcurrent); + } else { + fprintf(fp,"Atoms. Timestep: " BIGINT_FORMAT "\n",update->ntimestep); + } writer.write(header.c_str(), header.length()); } } diff --git a/src/dump_xyz.cpp b/src/dump_xyz.cpp index db929f9d6e..26e84d145d 100644 --- a/src/dump_xyz.cpp +++ b/src/dump_xyz.cpp @@ -131,11 +131,14 @@ void DumpXYZ::write_header(bigint n) { if (me == 0) { fprintf(fp,BIGINT_FORMAT "\n",n); - if (time_flag) + if (time_flag) { + double tcurrent = update->atime + + (update->ntimestep-update->atimestep) + update->dt; fprintf(fp,"Atoms. Timestep: " BIGINT_FORMAT " Time: %f\n", - update->ntimestep, update->atime); - else + update->ntimestep, tcurrent); + } else { fprintf(fp,"Atoms. Timestep: " BIGINT_FORMAT "\n",update->ntimestep); + } } } From cc4d7215f1778057c8cee1afa4f64963a76d7b93 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 21 Dec 2021 14:37:34 -0500 Subject: [PATCH 10/13] simplify. only output absolute time during MD. --- src/dump_xyz.cpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/dump_xyz.cpp b/src/dump_xyz.cpp index db929f9d6e..cb0a73eb81 100644 --- a/src/dump_xyz.cpp +++ b/src/dump_xyz.cpp @@ -130,12 +130,10 @@ int DumpXYZ::modify_param(int narg, char **arg) void DumpXYZ::write_header(bigint n) { if (me == 0) { - fprintf(fp,BIGINT_FORMAT "\n",n); - if (time_flag) - fprintf(fp,"Atoms. Timestep: " BIGINT_FORMAT " Time: %f\n", - update->ntimestep, update->atime); + if (time_flag && (update->whichflag == 1)) + fmt::print(fp,"{}\n Atoms. Timestep: {} Time: {:.6f}\n", update->ntimestep, update->atime); else - fprintf(fp,"Atoms. Timestep: " BIGINT_FORMAT "\n",update->ntimestep); + fmt::print(fp,"{}\n Atoms. Timestep: {}\n", update->ntimestep); } } From 6a442e1df43f1f029ae0387b6b4e3330e62f6719 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Tue, 21 Dec 2021 14:05:16 -0700 Subject: [PATCH 11/13] use compute_time() func in xyz output --- src/COMPRESS/dump_xyz_gz.cpp | 9 +++------ src/COMPRESS/dump_xyz_zstd.cpp | 9 +++------ src/dump.cpp | 1 + src/dump_xyz.cpp | 10 ++++------ 4 files changed, 11 insertions(+), 18 deletions(-) diff --git a/src/COMPRESS/dump_xyz_gz.cpp b/src/COMPRESS/dump_xyz_gz.cpp index a43a27ca7d..f017cb3b3e 100644 --- a/src/COMPRESS/dump_xyz_gz.cpp +++ b/src/COMPRESS/dump_xyz_gz.cpp @@ -94,14 +94,11 @@ void DumpXYZGZ::write_header(bigint ndump) { if (me == 0) { std::string header = fmt::format("{}\n", ndump); - if (time_flag) { - double tcurrent = update->atime + - (update->ntimestep-update->atimestep) + update->dt; + if (time_flag) fprintf(fp,"Atoms. Timestep: " BIGINT_FORMAT " Time: %f\n", - update->ntimestep, tcurrent); - } else { + update->ntimestep, compute_time()); + else header += fmt::format("Atoms. Timestep: {}\n", update->ntimestep); - } writer.write(header.c_str(), header.length()); } } diff --git a/src/COMPRESS/dump_xyz_zstd.cpp b/src/COMPRESS/dump_xyz_zstd.cpp index e955a0e3ba..46833f3c74 100644 --- a/src/COMPRESS/dump_xyz_zstd.cpp +++ b/src/COMPRESS/dump_xyz_zstd.cpp @@ -102,14 +102,11 @@ void DumpXYZZstd::write_header(bigint ndump) { if (me == 0) { std::string header = fmt::format("{}\n", ndump); - if (time_flag) { - double tcurrent = update->atime + - (update->ntimestep-update->atimestep) + update->dt; + if (time_flag) fprintf(fp,"Atoms. Timestep: " BIGINT_FORMAT " Time: %f\n", - update->ntimestep, tcurrent); - } else { + update->ntimestep, compute_time()); + else fprintf(fp,"Atoms. Timestep: " BIGINT_FORMAT "\n",update->ntimestep); - } writer.write(header.c_str(), header.length()); } } diff --git a/src/dump.cpp b/src/dump.cpp index df39f3738d..0f338e99b5 100644 --- a/src/dump.cpp +++ b/src/dump.cpp @@ -1132,6 +1132,7 @@ double Dump::compute_time() { return update->atime + (update->ntimestep - update->atimestep)*update->dt; } + /* ---------------------------------------------------------------------- return # of bytes of allocated memory ------------------------------------------------------------------------- */ diff --git a/src/dump_xyz.cpp b/src/dump_xyz.cpp index 26e84d145d..9943dd57a6 100644 --- a/src/dump_xyz.cpp +++ b/src/dump_xyz.cpp @@ -131,14 +131,12 @@ void DumpXYZ::write_header(bigint n) { if (me == 0) { fprintf(fp,BIGINT_FORMAT "\n",n); - if (time_flag) { - double tcurrent = update->atime + - (update->ntimestep-update->atimestep) + update->dt; + if (time_flag) fprintf(fp,"Atoms. Timestep: " BIGINT_FORMAT " Time: %f\n", - update->ntimestep, tcurrent); - } else { + update->ntimestep, compute_time()); + else fprintf(fp,"Atoms. Timestep: " BIGINT_FORMAT "\n",update->ntimestep); - } + } } } From d694b7cc1c3380c166f8a1da60520bf9152e0e8b Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 23 Dec 2021 14:34:49 -0500 Subject: [PATCH 12/13] recover compilation --- src/dump_xyz.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/dump_xyz.cpp b/src/dump_xyz.cpp index 9943dd57a6..949ddf1a11 100644 --- a/src/dump_xyz.cpp +++ b/src/dump_xyz.cpp @@ -137,7 +137,6 @@ void DumpXYZ::write_header(bigint n) else fprintf(fp,"Atoms. Timestep: " BIGINT_FORMAT "\n",update->ntimestep); } - } } /* ---------------------------------------------------------------------- */ From a653ee6b2ce37920c64596c03ed141d826933aa7 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 23 Dec 2021 15:19:17 -0500 Subject: [PATCH 13/13] recover failing unit tests and whitespace fixes --- src/EXTRA-DUMP/dump_xtc.cpp | 2 +- src/integrate.cpp | 4 +-- src/output.cpp | 40 +++++++++++++-------------- unittest/formats/test_dump_atom.cpp | 9 ++++++ unittest/formats/test_dump_custom.cpp | 10 +++++++ 5 files changed, 42 insertions(+), 23 deletions(-) diff --git a/src/EXTRA-DUMP/dump_xtc.cpp b/src/EXTRA-DUMP/dump_xtc.cpp index 94271f31a6..41b78ab64c 100644 --- a/src/EXTRA-DUMP/dump_xtc.cpp +++ b/src/EXTRA-DUMP/dump_xtc.cpp @@ -128,7 +128,7 @@ void DumpXTC::init_style() int idump; for (idump = 0; idump < output->ndump; idump++) if (strcmp(id,output->dump[idump]->id) == 0) break; - + if (output->mode_dump[idump] == 1) error->all(FLERR,"Cannot use every/time setting for dump xtc"); diff --git a/src/integrate.cpp b/src/integrate.cpp index c111532ee2..256291ed3b 100644 --- a/src/integrate.cpp +++ b/src/integrate.cpp @@ -120,7 +120,7 @@ void Integrate::ev_setup() (1) computes that need energy/virial info on this timestep (2) time dumps that may need per-atom compute info on this timestep NOTE: inefficient to add all per-atom eng/virial computes - but don't know which ones the dump needs + but don't know which ones the dump needs see NOTE in output.cpp invoke matchstep() on all timestep-dependent computes to clear their arrays eflag: set any or no bits @@ -140,7 +140,7 @@ void Integrate::ev_set(bigint ntimestep) int i,flag; int tdflag = 0; - if (output->any_time_dumps && + if (output->any_time_dumps && output->next_time_dump_any == ntimestep) tdflag = 1; flag = 0; diff --git a/src/output.cpp b/src/output.cpp index d828b5a417..b7b1a0d8c1 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -199,9 +199,9 @@ void Output::setup(int memflag) // (2) and (3) only apply for non-variable dump intervals // finally, do not write if same snapshot written previously, // i.e. on last timestep of previous run - + int writeflag = 0; - + if (last_dump[idump] < 0 && dump[idump]->first_flag == 1) writeflag = 1; if (mode_dump[idump] == 0) { @@ -209,7 +209,7 @@ void Output::setup(int memflag) writeflag = 1; } else { if (every_time_dump[idump] > 0.0) { - double tcurrent = update->atime + + double tcurrent = update->atime + (ntimestep - update->atimestep) * update->dt; double remainder = fmod(tcurrent,every_time_dump[idump]); if ((remainder < EPSDT*update->dt) || @@ -234,7 +234,7 @@ void Output::setup(int memflag) if (writeflag || last_dump[idump] < 0) calculate_next_dump(0,idump,ntimestep); - // if dump not written now, use addstep_compute_all() + // if dump not written now, use addstep_compute_all() // since don't know what computes the dump will invoke if (mode_dump[idump] == 0 && @@ -354,7 +354,7 @@ void Output::setup(int memflag) if (next_dump[idump] == ntimestep) { if (last_dump[idump] == ntimestep) continue; - if (mode_dump[idump] == 0 && + if (mode_dump[idump] == 0 && (dump[idump]->clearstep || var_dump[idump])) modify->clearstep_compute(); @@ -365,7 +365,7 @@ void Output::setup(int memflag) last_dump[idump] = ntimestep; calculate_next_dump(1,idump,ntimestep); - if (mode_dump[idump] == 0 && + if (mode_dump[idump] == 0 && (dump[idump]->clearstep || var_dump[idump])) modify->addstep_compute(next_dump[idump]); } @@ -490,7 +490,7 @@ void Output::setup(int memflag) if (which == 0) next_dump[idump] = (ntimestep/every_dump[idump])*every_dump[idump] + every_dump[idump]; - else if (which == 1) + else if (which == 1) next_dump[idump] += every_dump[idump]; } else { @@ -507,7 +507,7 @@ void Output::setup(int memflag) bigint nextdump; double nexttime; - double tcurrent = update->atime + + double tcurrent = update->atime + (ntimestep - update->atimestep) * update->dt; if (every_time_dump[idump] > 0.0) { @@ -517,27 +517,27 @@ void Output::setup(int memflag) // which = 2: no change to previous nexttime (only timestep has changed) if (which == 0) - nexttime = static_cast (tcurrent/every_time_dump[idump]) * + nexttime = static_cast (tcurrent/every_time_dump[idump]) * every_time_dump[idump] + every_time_dump[idump]; - else if (which == 1) + else if (which == 1) nexttime = next_time_dump[idump] + every_time_dump[idump]; else if (which == 2) nexttime = next_time_dump[idump]; - nextdump = ntimestep + - static_cast ((nexttime - tcurrent - EPSDT*update->dt) / + nextdump = ntimestep + + static_cast ((nexttime - tcurrent - EPSDT*update->dt) / update->dt) + 1; // if delta is too small to reach next timestep, use multiple of delta if (nextdump == ntimestep) { - double tnext = update->atime + + double tnext = update->atime + (ntimestep+1 - update->atimestep) * update->dt; - int multiple = static_cast + int multiple = static_cast ((tnext - nexttime) / every_time_dump[idump]); nexttime = nexttime + (multiple+1)*every_time_dump[idump]; - nextdump = ntimestep + - static_cast ((nexttime - tcurrent - EPSDT*update->dt) / + nextdump = ntimestep + + static_cast ((nexttime - tcurrent - EPSDT*update->dt) / update->dt) + 1; } @@ -548,14 +548,14 @@ void Output::setup(int memflag) if (which < 2 || next_time_dump[idump] < 0.0) { nexttime = input->variable->compute_equal(ivar_dump[idump]); - } else + } else nexttime = next_time_dump[idump]; if (nexttime <= tcurrent) error->all(FLERR,"Dump every/time variable returned a bad time"); - nextdump = ntimestep + - static_cast ((nexttime - tcurrent - EPSDT*update->dt) / + nextdump = ntimestep + + static_cast ((nexttime - tcurrent - EPSDT*update->dt) / update->dt) + 1; if (nextdump <= ntimestep) error->all(FLERR,"Dump every/time variable too small for next timestep"); @@ -704,7 +704,7 @@ void Output::reset_dt() // reset next_dump but do not change next_time_dump, 2 arg for reset_dt() // do not invoke for a dump already scheduled for this step // since timestep change affects next step - + if (next_dump[idump] != ntimestep) calculate_next_dump(2,idump,update->ntimestep); diff --git a/unittest/formats/test_dump_atom.cpp b/unittest/formats/test_dump_atom.cpp index 1d00e00610..a73204fb92 100644 --- a/unittest/formats/test_dump_atom.cpp +++ b/unittest/formats/test_dump_atom.cpp @@ -74,6 +74,13 @@ public: END_HIDE_OUTPUT(); } + void close_dump() + { + BEGIN_HIDE_OUTPUT(); + command("undump id"); + END_HIDE_OUTPUT(); + } + void generate_text_and_binary_dump(std::string text_file, std::string binary_file, std::string dump_modify_options, int ntimesteps) { @@ -505,6 +512,7 @@ TEST_F(DumpAtomTest, rerun) ASSERT_FILE_EXISTS(dump_file); ASSERT_EQ(count_lines(dump_file), 82); continue_dump(1); + close_dump(); lmp->output->thermo->evaluate_keyword("pe", &pe_2); ASSERT_FILE_EXISTS(dump_file); ASSERT_EQ(count_lines(dump_file), 123); @@ -532,6 +540,7 @@ TEST_F(DumpAtomTest, rerun_bin) lmp->output->thermo->evaluate_keyword("pe", &pe_1); ASSERT_FILE_EXISTS(dump_file); continue_dump(1); + close_dump(); lmp->output->thermo->evaluate_keyword("pe", &pe_2); ASSERT_FILE_EXISTS(dump_file); HIDE_OUTPUT([&] { diff --git a/unittest/formats/test_dump_custom.cpp b/unittest/formats/test_dump_custom.cpp index 921b584217..434acf462c 100644 --- a/unittest/formats/test_dump_custom.cpp +++ b/unittest/formats/test_dump_custom.cpp @@ -73,6 +73,13 @@ public: END_HIDE_OUTPUT(); } + void close_dump() + { + BEGIN_HIDE_OUTPUT(); + command("undump id"); + END_HIDE_OUTPUT(); + } + void generate_text_and_binary_dump(std::string text_file, std::string binary_file, std::string fields, std::string dump_modify_options, int ntimesteps) @@ -330,6 +337,7 @@ TEST_F(DumpCustomTest, rerun) ASSERT_FILE_EXISTS(dump_file); ASSERT_EQ(count_lines(dump_file), 82); continue_dump(1); + close_dump(); lmp->output->thermo->evaluate_keyword("pe", &pe_2); ASSERT_FILE_EXISTS(dump_file); ASSERT_EQ(count_lines(dump_file), 123); @@ -338,6 +346,7 @@ TEST_F(DumpCustomTest, rerun) }); lmp->output->thermo->evaluate_keyword("pe", &pe_rerun); ASSERT_DOUBLE_EQ(pe_1, pe_rerun); + HIDE_OUTPUT([&] { command(fmt::format("rerun {} first 2 last 2 every 1 post yes dump x y z", dump_file)); }); @@ -359,6 +368,7 @@ TEST_F(DumpCustomTest, rerun_bin) lmp->output->thermo->evaluate_keyword("pe", &pe_1); ASSERT_FILE_EXISTS(dump_file); continue_dump(1); + close_dump(); lmp->output->thermo->evaluate_keyword("pe", &pe_2); ASSERT_FILE_EXISTS(dump_file); HIDE_OUTPUT([&] {