diff --git a/src/compute.h b/src/compute.h index 90b88c7d1d..8151bd8ee0 100644 --- a/src/compute.h +++ b/src/compute.h @@ -30,6 +30,7 @@ class Compute : protected Pointers { INVOKED_ARRAY = 1<<2, INVOKED_PERATOM = 1<<3, INVOKED_LOCAL = 1<<4, + INVOKED_PERGRID = 1<<5, }; // clang-format on static int instance_total; // # of Compute classes ever instantiated @@ -61,6 +62,8 @@ class Compute : protected Pointers { int size_local_rows; // rows in local vector or array int size_local_cols; // 0 = vector, N = columns in local array + int pergrid_flag; // 0/1 if compute_pergrid() function exists + int extscalar; // 0/1 if global scalar is intensive/extensive int extvector; // 0/1/-1 if global vector is all int/ext/extlist int *extlist; // list of 0/1 int/ext for each vec component @@ -91,6 +94,7 @@ class Compute : protected Pointers { bigint invoked_array; // ditto for compute_array() bigint invoked_peratom; // ditto for compute_peratom() bigint invoked_local; // ditto for compute_local() + bigint invoked_pergrid; // ditto for compute_grid() double dof; // degrees-of-freedom for temperature @@ -118,6 +122,7 @@ class Compute : protected Pointers { virtual void compute_array() {} virtual void compute_peratom() {} virtual void compute_local() {} + virtual void compute_pergrid() {} virtual void set_arrays(int) {} virtual int pack_forward_comm(int, int *, double *, int, int *) { return 0; } diff --git a/src/dump_grid.cpp b/src/dump_grid.cpp new file mode 100644 index 0000000000..2f6c510611 --- /dev/null +++ b/src/dump_grid.cpp @@ -0,0 +1,863 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "dump_grid.h" + +#include "arg_info.h" +#include "compute.h" +#include "domain.h" +#include "error.h" +#include "fix.h" +#include "memory.h" +#include "modify.h" +#include "region.h" +#include "update.h" + +#include + +using namespace LAMMPS_NS; + +// customize by adding keyword +// also customize compute_atom_property.cpp + +enum{COMPUTE,FIX}; + +#define ONEFIELD 32 +#define DELTA 1048576 + +/* ---------------------------------------------------------------------- */ + +DumpGrid::DumpGrid(LAMMPS *lmp, int narg, char **arg) : + Dump(lmp, narg, arg) +{ + if (narg == 5) error->all(FLERR,"No dump grid arguments specified"); + + clearstep = 1; + + nevery = utils::inumeric(FLERR,arg[3],false,lmp); + if (nevery <= 0) error->all(FLERR,"Illegal dump grid command"); + + // expand args if any have wildcard character "*" + // ok to include trailing optional args, + // so long as they do not have "*" between square brackets + // nfield may be shrunk below if extra optional args exist + + expand = 0; + nfield = nargnew = utils::expand_args(FLERR,narg-5,&arg[5],1,earg,lmp); + if (earg != &arg[5]) expand = 1; + + // allocate field vectors + + pack_choice = new FnPtrPack[nfield]; + vtype = new int[nfield]; + memory->create(field2index,nfield,"dump:field2index"); + memory->create(argindex,nfield,"dump:argindex"); + + buffer_allow = 1; + buffer_flag = 1; + + // computes and fixes which the dump accesses + + ncompute = 0; + nfix = 0; + + // process attributes + // ioptional = start of additional optional args in expanded args + + ioptional = parse_fields(nfield,earg); + + if (ioptional < nfield && + strcmp(style,"image") != 0 && strcmp(style,"movie") != 0) + error->all(FLERR,"Invalid attribute {} in dump {} command",earg[ioptional],style); + + // noptional = # of optional args + // reset nfield to subtract off optional args + // reset ioptional to what it would be in original arg list + // only dump image and dump movie styles process optional args, + // they do not use expanded earg list + + int noptional = nfield - ioptional; + nfield -= noptional; + size_one = nfield; + ioptional = narg - noptional; + + // setup format strings + + vformat = new char*[nfield]; + std::string cols; + + cols.clear(); + for (int i = 0; i < nfield; i++) { + if (vtype[i] == Dump::INT) cols += "%d "; + else if (vtype[i] == Dump::DOUBLE) cols += "%g "; + else if (vtype[i] == Dump::STRING) cols += "%s "; + else if (vtype[i] == Dump::BIGINT) cols += BIGINT_FORMAT " "; + vformat[i] = nullptr; + } + cols.resize(cols.size()-1); + format_default = utils::strdup(cols); + + format_column_user = new char*[nfield]; + for (int i = 0; i < nfield; i++) format_column_user[i] = nullptr; + + // setup column string + + cols.clear(); + keyword_user.resize(nfield); + for (int iarg = 0; iarg < nfield; iarg++) { + key2col[earg[iarg]] = iarg; + keyword_user[iarg].clear(); + if (cols.size()) cols += " "; + cols += earg[iarg]; + } + columns_default = utils::strdup(cols); +} + +/* ---------------------------------------------------------------------- */ + +DumpGrid::~DumpGrid() +{ + // if wildcard expansion occurred, free earg memory from expand_args() + // could not do in constructor, b/c some derived classes process earg + + if (expand) { + for (int i = 0; i < nargnew; i++) delete[] earg[i]; + memory->sfree(earg); + } + + delete[] pack_choice; + delete[] vtype; + memory->destroy(field2index); + memory->destroy(argindex); + + delete[] idregion; + + for (int i = 0; i < ncompute; i++) delete[] id_compute[i]; + memory->sfree(id_compute); + delete[] compute; + + for (int i = 0; i < nfix; i++) delete[] id_fix[i]; + memory->sfree(id_fix); + delete[] fix; + + if (vformat) { + for (int i = 0; i < nfield; i++) delete[] vformat[i]; + delete[] vformat; + } + + if (format_column_user) { + for (int i = 0; i < nfield; i++) delete[] format_column_user[i]; + delete[] format_column_user; + } + + delete[] columns_default; + delete[] columns; +} + +/* ---------------------------------------------------------------------- */ + +void DumpGrid::init_style() +{ + // assemble ITEMS: column string from defaults and user values + + delete[] columns; + std::string combined; + int icol = 0; + for (auto item : utils::split_words(columns_default)) { + if (combined.size()) combined += " "; + if (keyword_user[icol].size()) combined += keyword_user[icol]; + else combined += item; + ++icol; + } + columns = utils::strdup(combined); + + // format = copy of default or user-specified line format + + delete[] format; + if (format_line_user) format = utils::strdup(format_line_user); + else format = utils::strdup(format_default); + + // tokenize the format string and add space at end of each format element + // if user-specified int/float format exists, use it instead + // if user-specified column format exists, use it instead + // lo priority = line, medium priority = int/float, hi priority = column + + auto words = utils::split_words(format); + if ((int) words.size() < nfield) + error->all(FLERR,"Dump_modify format line is too short"); + + int i=0; + for (const auto &word : words) { + delete[] vformat[i]; + + if (format_column_user[i]) + vformat[i] = utils::strdup(std::string(format_column_user[i]) + " "); + else if (vtype[i] == Dump::INT && format_int_user) + vformat[i] = utils::strdup(std::string(format_int_user) + " "); + else if (vtype[i] == Dump::DOUBLE && format_float_user) + vformat[i] = utils::strdup(std::string(format_float_user) + " "); + else if (vtype[i] == Dump::BIGINT && format_bigint_user) + vformat[i] = utils::strdup(std::string(format_bigint_user) + " "); + else vformat[i] = utils::strdup(word + " "); + + // remove trailing blank on last column's format + if (i == nfield-1) vformat[i][strlen(vformat[i])-1] = '\0'; + + ++i; + } + + // setup boundary string + + domain->boundary_string(boundstr); + + // setup function ptrs + + if (binary && domain->triclinic == 0) + header_choice = &DumpGrid::header_binary; + else if (binary && domain->triclinic == 1) + header_choice = &DumpGrid::header_binary_triclinic; + else if (!binary && domain->triclinic == 0) + header_choice = &DumpGrid::header_item; + else if (!binary && domain->triclinic == 1) + header_choice = &DumpGrid::header_item_triclinic; + + if (binary) write_choice = &DumpGrid::write_binary; + else if (buffer_flag == 1) write_choice = &DumpGrid::write_string; + else write_choice = &DumpGrid::write_lines; + + // find current ptr for each compute and fix + // check that fix frequency is acceptable + + for (i = 0; i < ncompute; i++) { + compute[i] = modify->get_compute_by_id(id_compute[i]); + if (!compute[i]) error->all(FLERR,"Could not find dump grid compute ID {}",id_compute[i]); + } + + for (i = 0; i < nfix; i++) { + fix[i] = modify->get_fix_by_id(id_fix[i]); + if (!fix[i]) error->all(FLERR,"Could not find dump grid fix ID {}", id_fix[i]); + if (nevery % fix[i]->peratom_freq) + error->all(FLERR,"Dump grid and fix not computed at compatible times"); + } + + // check validity of region + + if (idregion && !domain->get_region_by_id(idregion)) + error->all(FLERR,"Region {} for dump grid does not exist", idregion); + + // open single file, one time only + + if (multifile == 0) openfile(); +} + +/* ---------------------------------------------------------------------- */ + +void DumpGrid::write_header(bigint ndump) +{ + if (multiproc) (this->*header_choice)(ndump); + else if (me == 0) (this->*header_choice)(ndump); +} + +/* ---------------------------------------------------------------------- */ + +void DumpGrid::format_magic_string_binary() +{ + // use negative ntimestep as marker for new format + bigint fmtlen = strlen(MAGIC_STRING); + bigint marker = -fmtlen; + fwrite(&marker, sizeof(bigint), 1, fp); + fwrite(MAGIC_STRING, sizeof(char), fmtlen, fp); +} + +/* ---------------------------------------------------------------------- */ + +void DumpGrid::format_endian_binary() +{ + int endian = ENDIAN; + fwrite(&endian, sizeof(int), 1, fp); +} + +/* ---------------------------------------------------------------------- */ + +void DumpGrid::format_revision_binary() +{ + int revision = FORMAT_REVISION; + fwrite(&revision, sizeof(int), 1, fp); +} + +/* ---------------------------------------------------------------------- */ + +void DumpGrid::header_unit_style_binary() +{ + int len = 0; + if (unit_flag && !unit_count) { + ++unit_count; + len = strlen(update->unit_style); + fwrite(&len, sizeof(int), 1, fp); + fwrite(update->unit_style, sizeof(char), len, fp); + } else { + fwrite(&len, sizeof(int), 1, fp); + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpGrid::header_columns_binary() +{ + int len = strlen(columns); + fwrite(&len, sizeof(int), 1, fp); + fwrite(columns, sizeof(char), len, fp); +} + +/* ---------------------------------------------------------------------- */ + +void DumpGrid::header_time_binary() +{ + char flag = time_flag ? 1 : 0; + fwrite(&flag, sizeof(char), 1, fp); + + if (time_flag) { + double t = compute_time(); + fwrite(&t, sizeof(double), 1, fp); + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpGrid::header_format_binary() +{ + format_magic_string_binary(); + format_endian_binary(); + format_revision_binary(); +} + +/* ---------------------------------------------------------------------- */ + +void DumpGrid::header_binary(bigint ndump) +{ + header_format_binary(); + + fwrite(&update->ntimestep,sizeof(bigint),1,fp); + fwrite(&ndump,sizeof(bigint),1,fp); + fwrite(&domain->triclinic,sizeof(int),1,fp); + fwrite(&domain->boundary[0][0],6*sizeof(int),1,fp); + fwrite(&boxxlo,sizeof(double),1,fp); + fwrite(&boxxhi,sizeof(double),1,fp); + fwrite(&boxylo,sizeof(double),1,fp); + fwrite(&boxyhi,sizeof(double),1,fp); + fwrite(&boxzlo,sizeof(double),1,fp); + fwrite(&boxzhi,sizeof(double),1,fp); + fwrite(&nfield,sizeof(int),1,fp); + + header_unit_style_binary(); + header_time_binary(); + header_columns_binary(); + + if (multiproc) fwrite(&nclusterprocs,sizeof(int),1,fp); + else fwrite(&nprocs,sizeof(int),1,fp); +} + +/* ---------------------------------------------------------------------- */ + +void DumpGrid::header_binary_triclinic(bigint ndump) +{ + header_format_binary(); + + fwrite(&update->ntimestep,sizeof(bigint),1,fp); + fwrite(&ndump,sizeof(bigint),1,fp); + fwrite(&domain->triclinic,sizeof(int),1,fp); + fwrite(&domain->boundary[0][0],6*sizeof(int),1,fp); + fwrite(&boxxlo,sizeof(double),1,fp); + fwrite(&boxxhi,sizeof(double),1,fp); + fwrite(&boxylo,sizeof(double),1,fp); + fwrite(&boxyhi,sizeof(double),1,fp); + fwrite(&boxzlo,sizeof(double),1,fp); + fwrite(&boxzhi,sizeof(double),1,fp); + fwrite(&boxxy,sizeof(double),1,fp); + fwrite(&boxxz,sizeof(double),1,fp); + fwrite(&boxyz,sizeof(double),1,fp); + fwrite(&nfield,sizeof(int),1,fp); + + header_unit_style_binary(); + header_time_binary(); + header_columns_binary(); + + if (multiproc) fwrite(&nclusterprocs,sizeof(int),1,fp); + else fwrite(&nprocs,sizeof(int),1,fp); +} + +/* ---------------------------------------------------------------------- */ + +void DumpGrid::header_item(bigint ndump) +{ + if (unit_flag && !unit_count) { + ++unit_count; + fmt::print(fp,"ITEM: UNITS\n{}\n",update->unit_style); + } + if (time_flag) fmt::print(fp,"ITEM: TIME\n{:.16}\n",compute_time()); + + fmt::print(fp,"ITEM: TIMESTEP\n{}\n" + "ITEM: NUMBER OF ATOMS\n{}\n", + update->ntimestep, ndump); + + fmt::print(fp,"ITEM: BOX BOUNDS {}\n" + "{:>1.16e} {:>1.16e}\n" + "{:>1.16e} {:>1.16e}\n" + "{:>1.16e} {:>1.16e}\n", + boundstr,boxxlo,boxxhi,boxylo,boxyhi,boxzlo,boxzhi); + + fmt::print(fp,"ITEM: ATOMS {}\n",columns); +} + +/* ---------------------------------------------------------------------- */ + +void DumpGrid::header_item_triclinic(bigint ndump) +{ + if (unit_flag && !unit_count) { + ++unit_count; + fmt::print(fp,"ITEM: UNITS\n{}\n",update->unit_style); + } + if (time_flag) fmt::print(fp,"ITEM: TIME\n{:.16}\n",compute_time()); + + fmt::print(fp,"ITEM: TIMESTEP\n{}\n" + "ITEM: NUMBER OF ATOMS\n{}\n", + update->ntimestep, ndump); + + fmt::print(fp,"ITEM: BOX BOUNDS xy xz yz {}\n" + "{:>1.16e} {:>1.16e} {:>1.16e}\n" + "{:>1.16e} {:>1.16e} {:>1.16e}\n" + "{:>1.16e} {:>1.16e} {:>1.16e}\n", + boundstr,boxxlo,boxxhi,boxxy,boxylo,boxyhi,boxxz,boxzlo,boxzhi,boxyz); + + fmt::print(fp,"ITEM: ATOMS {}\n",columns); +} + +/* ---------------------------------------------------------------------- */ + +int DumpGrid::count() +{ + int i; + + // grow choose arrays if needed + // NOTE: needs to change + + /* + const int nlocal = atom->nlocal; + if (atom->nmax > maxlocal) { + maxlocal = atom->nmax; + + memory->destroy(choose); + memory->destroy(dchoose); + memory->destroy(clist); + memory->create(choose,maxlocal,"dump:choose"); + memory->create(dchoose,maxlocal,"dump:dchoose"); + memory->create(clist,maxlocal,"dump:clist"); + } + */ + + // invoke Computes for per-grid quantities + // only if within a run or minimize + // else require that computes are current + // this prevents a compute from being invoked by the WriteDump class + + if (ncompute) { + if (update->whichflag == 0) { + for (i = 0; i < ncompute; i++) + if (compute[i]->invoked_pergrid != update->ntimestep) + error->all(FLERR,"Compute used in dump between runs is not current"); + } else { + for (i = 0; i < ncompute; i++) { + if (!(compute[i]->invoked_flag & Compute::INVOKED_PERGRID)) { + compute[i]->compute_pergrid(); + compute[i]->invoked_flag |= Compute::INVOKED_PERGRID; + } + } + } + } + + // choose all local grid pts for output + // NOTE: this needs to change + + //for (i = 0; i < nlocal; i++) choose[i] = 1; + + // un-choose if not in region + // NOTE: this needs to change + + if (idregion) { + auto region = domain->get_region_by_id(idregion); + region->prematch(); + /* + double **x = atom->x; + for (i = 0; i < nlocal; i++) + if (choose[i] && region->match(x[i][0],x[i][1],x[i][2]) == 0) + choose[i] = 0; + */ + } + + // compress choose flags into clist + // nchoose = # of selected atoms + // clist[i] = local index of each selected atom + // NOTE: this neds to change + + nchoose = 0; + /* + for (i = 0; i < nlocal; i++) + if (choose[i]) clist[nchoose++] = i; + */ + + return nchoose; +} + +/* ---------------------------------------------------------------------- */ + +void DumpGrid::pack(tagint *ids) +{ + for (int n = 0; n < size_one; n++) (this->*pack_choice[n])(n); + // NOTE: this needs to be grid IDs ? + /* + if (ids) { + tagint *tag = atom->tag; + for (int i = 0; i < nchoose; i++) + ids[i] = tag[clist[i]]; + } + */ +} + +/* ---------------------------------------------------------------------- + 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 +------------------------------------------------------------------------- */ + +int DumpGrid::convert_string(int n, double *mybuf) +{ + int i,j; + + int offset = 0; + int m = 0; + for (i = 0; i < n; i++) { + if (offset + nfield*ONEFIELD > maxsbuf) { + if ((bigint) maxsbuf + DELTA > MAXSMALLINT) return -1; + maxsbuf += DELTA; + memory->grow(sbuf,maxsbuf,"dump:sbuf"); + } + + for (j = 0; j < nfield; j++) { + if (vtype[j] == Dump::INT) + offset += sprintf(&sbuf[offset],vformat[j],static_cast (mybuf[m])); + else if (vtype[j] == Dump::DOUBLE) + offset += sprintf(&sbuf[offset],vformat[j],mybuf[m]); + else if (vtype[j] == Dump::BIGINT) + offset += sprintf(&sbuf[offset],vformat[j], + static_cast (mybuf[m])); + m++; + } + offset += sprintf(&sbuf[offset],"\n"); + } + + return offset; +} + +/* ---------------------------------------------------------------------- */ + +void DumpGrid::write_data(int n, double *mybuf) +{ + (this->*write_choice)(n,mybuf); +} + +/* ---------------------------------------------------------------------- */ + +void DumpGrid::write_binary(int n, double *mybuf) +{ + n *= size_one; + fwrite(&n,sizeof(int),1,fp); + fwrite(mybuf,sizeof(double),n,fp); +} + +/* ---------------------------------------------------------------------- */ + +void DumpGrid::write_string(int n, double *mybuf) +{ + if (mybuf) + fwrite(mybuf,sizeof(char),n,fp); +} + +/* ---------------------------------------------------------------------- */ + +void DumpGrid::write_lines(int n, double *mybuf) +{ + int i,j; + + int m = 0; + for (i = 0; i < n; i++) { + for (j = 0; j < nfield; j++) { + if (vtype[j] == Dump::INT) fprintf(fp,vformat[j],static_cast (mybuf[m])); + else if (vtype[j] == Dump::DOUBLE) fprintf(fp,vformat[j],mybuf[m]); + else if (vtype[j] == Dump::BIGINT) + fprintf(fp,vformat[j],static_cast (mybuf[m])); + m++; + } + fprintf(fp,"\n"); + } +} + +/* ---------------------------------------------------------------------- */ + +int DumpGrid::parse_fields(int narg, char **arg) +{ + // customize by adding to if statement + + for (int iarg = 0; iarg < narg; iarg++) { + + int n,flag,cols; + ArgInfo argi(arg[iarg], ArgInfo::COMPUTE | ArgInfo::FIX); + argindex[iarg] = argi.get_index1(); + auto name = argi.get_name(); + Compute *icompute = nullptr; + Fix *ifix = nullptr; + + switch (argi.get_type()) { + + case ArgInfo::UNKNOWN: + error->all(FLERR,"Invalid attribute in dump grid command"); + break; + + // compute value = c_ID + // if no trailing [], then arg is set to 0, else arg is int between [] + + case ArgInfo::COMPUTE: + pack_choice[iarg] = &DumpGrid::pack_compute; + vtype[iarg] = Dump::DOUBLE; + + icompute = modify->get_compute_by_id(name); + if (!icompute) error->all(FLERR,"Could not find dump grid compute ID: {}",name); + if (icompute->pergrid_flag == 0) + error->all(FLERR,"Dump grid compute {} does not compute per-grid info",name); + /* + if (argi.get_dim() == 0 && icompute->size_pergrid_cols > 0) + error->all(FLERR,"Dump grid compute {} does not calculate per-grid vector",name); + if (argi.get_dim() > 0 && icompute->size_pergrid_cols == 0) + error->all(FLERR,"Dump grid compute {} does not calculate per-grid array",name); + if (argi.get_dim() > 0 && argi.get_index1() > icompute->size_pergrid_cols) + error->all(FLERR,"Dump grid compute {} vector is accessed out-of-range",name); + */ + + field2index[iarg] = add_compute(name); + break; + + // fix value = f_ID + // if no trailing [], then arg is set to 0, else arg is between [] + + case ArgInfo::FIX: + pack_choice[iarg] = &DumpGrid::pack_fix; + vtype[iarg] = Dump::DOUBLE; + + ifix = modify->get_fix_by_id(name); + if (!ifix) error->all(FLERR,"Could not find dump grid fix ID: {}",name); + if (ifix->pergrid_flag == 0) + error->all(FLERR,"Dump grid fix {} does not compute per-atom info",name); + + /* + if (argi.get_dim() == 0 && ifix->size_pergrid_cols > 0) + error->all(FLERR,"Dump grid fix {} does not compute per-atom vector",name); + if (argi.get_dim() > 0 && ifix->size_pergrid_cols == 0) + error->all(FLERR,"Dump grid fix {} does not compute per-atom array",name); + if (argi.get_dim() > 0 && argi.get_index1() > ifix->size_pergrid_cols) + error->all(FLERR,"Dump grid fix {} vector is accessed out-of-range",name); + */ + + field2index[iarg] = add_fix(name); + break; + + // no match + + default: + return iarg; + break; + } + } + + return narg; +} + +/* ---------------------------------------------------------------------- + add Compute to list of Compute objects used by dump + return index of where this Compute is in list + if already in list, do not add, just return index, else add to list +------------------------------------------------------------------------- */ + +int DumpGrid::add_compute(const char *id) +{ + int icompute; + for (icompute = 0; icompute < ncompute; icompute++) + if (strcmp(id,id_compute[icompute]) == 0) break; + if (icompute < ncompute) return icompute; + + id_compute = (char **) + memory->srealloc(id_compute,(ncompute+1)*sizeof(char *),"dump:id_compute"); + delete[] compute; + compute = new Compute*[ncompute+1]; + + id_compute[ncompute] = utils::strdup(id); + ncompute++; + return ncompute-1; +} + +/* ---------------------------------------------------------------------- + add Fix to list of Fix objects used by dump + return index of where this Fix is in list + if already in list, do not add, just return index, else add to list +------------------------------------------------------------------------- */ + +int DumpGrid::add_fix(const char *id) +{ + int ifix; + for (ifix = 0; ifix < nfix; ifix++) + if (strcmp(id,id_fix[ifix]) == 0) break; + if (ifix < nfix) return ifix; + + id_fix = (char **) + memory->srealloc(id_fix,(nfix+1)*sizeof(char *),"dump:id_fix"); + delete[] fix; + fix = new Fix*[nfix+1]; + + id_fix[nfix] = utils::strdup(id); + nfix++; + return nfix-1; +} + +/* ---------------------------------------------------------------------- */ + +int DumpGrid::modify_param(int narg, char **arg) +{ + if (strcmp(arg[0],"region") == 0) { + if (narg < 2) error->all(FLERR,"Illegal dump_modify command"); + if (strcmp(arg[1],"none") == 0) { + delete[] idregion; + idregion = nullptr; + } else { + delete[] idregion; + if (!domain->get_region_by_id(arg[1])) + error->all(FLERR,"Dump_modify region {} does not exist", arg[1]); + idregion = utils::strdup(arg[1]); + } + return 2; + } + + if (strcmp(arg[0],"format") == 0) { + if (narg < 2) error->all(FLERR,"Illegal dump_modify command"); + + if (strcmp(arg[1],"none") == 0) { + // just clear format_column_user allocated by this dump child class + for (int i = 0; i < nfield; i++) { + delete[] format_column_user[i]; + format_column_user[i] = nullptr; + } + return 2; + } + + if (narg < 3) error->all(FLERR,"Illegal dump_modify command"); + + if (strcmp(arg[1],"int") == 0) { + delete[] format_int_user; + format_int_user = utils::strdup(arg[2]); + delete[] format_bigint_user; + int n = strlen(format_int_user) + 8; + format_bigint_user = new char[n]; + // replace "d" in format_int_user with bigint format specifier + // use of &str[1] removes leading '%' from BIGINT_FORMAT string + char *ptr = strchr(format_int_user,'d'); + if (ptr == nullptr) + error->all(FLERR,"Dump_modify int format does not contain d character"); + char str[8]; + sprintf(str,"%s",BIGINT_FORMAT); + *ptr = '\0'; + sprintf(format_bigint_user,"%s%s%s",format_int_user,&str[1],ptr+1); + *ptr = 'd'; + + } else if (strcmp(arg[1],"float") == 0) { + delete[] format_float_user; + format_float_user = utils::strdup(arg[2]); + + } else { + int i = utils::inumeric(FLERR,arg[1],false,lmp) - 1; + if (i < 0 || i >= nfield) + error->all(FLERR,"Illegal dump_modify command"); + delete[] format_column_user[i]; + format_column_user[i] = utils::strdup(arg[2]); + } + return 3; + } + + return 0; +} + +/* ---------------------------------------------------------------------- + return # of bytes of allocated memory in buf, choose, variable arrays +------------------------------------------------------------------------- */ + +double DumpGrid::memory_usage() +{ + double bytes = Dump::memory_usage(); + bytes += memory->usage(choose,maxlocal); + bytes += memory->usage(dchoose,maxlocal); + bytes += memory->usage(clist,maxlocal); + return bytes; +} + +/* ---------------------------------------------------------------------- + extraction of Compute and Fix data +------------------------------------------------------------------------- */ + +void DumpGrid::pack_compute(int n) +{ + double *vector = compute[field2index[n]]->vector_atom; + double **array = compute[field2index[n]]->array_atom; + int index = argindex[n]; + + if (index == 0) { + for (int i = 0; i < nchoose; i++) { + buf[n] = vector[clist[i]]; + n += size_one; + } + } else { + index--; + for (int i = 0; i < nchoose; i++) { + buf[n] = array[clist[i]][index]; + n += size_one; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpGrid::pack_fix(int n) +{ + double *vector = fix[field2index[n]]->vector_atom; + double **array = fix[field2index[n]]->array_atom; + int index = argindex[n]; + + if (index == 0) { + for (int i = 0; i < nchoose; i++) { + buf[n] = vector[clist[i]]; + n += size_one; + } + } else { + index--; + for (int i = 0; i < nchoose; i++) { + buf[n] = array[clist[i]][index]; + n += size_one; + } + } +} + diff --git a/src/dump_grid.h b/src/dump_grid.h new file mode 100644 index 0000000000..95714bd552 --- /dev/null +++ b/src/dump_grid.h @@ -0,0 +1,119 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef DUMP_CLASS +// clang-format off +DumpStyle(grid,DumpGrid); +// clang-format on +#else + +#ifndef LMP_DUMP_GRID_H +#define LMP_DUMP_GRID_H + +#include "dump.h" + +namespace LAMMPS_NS { + +class DumpGrid : public Dump { + public: + DumpGrid(class LAMMPS *, int, char **); + ~DumpGrid() override; + + const char *MAGIC_STRING = "DUMPGRID"; + const int FORMAT_REVISION = 0x0002; + const int ENDIAN = 0x0001; + + protected: + int nevery; // dump frequency for output + char *idregion; // region ID, nullptr if no region + + int expand; // flag for whether field args were expanded + char **earg; // field names with wildcard expansion + int nargnew; // size of earg + // + int *vtype; // type of each vector (INT, DOUBLE) + char **vformat; // format string for each vector element + // + char *columns; // column labels + char *columns_default; + + int nchoose; // # of selected atoms + int maxlocal; // size of atom selection and variable arrays + int *choose; // local indices of selected atoms + double *dchoose; // value for each atom to threshold against + int *clist; // compressed list of indices of selected atoms + + int nfield; // # of keywords listed by user + int ioptional; // index of start of optional args + // + int *field2index; // which compute,fix,variable,custom calcs this field + int *argindex; // index into compute,fix,custom per-atom data + // 0 for per-atom vector, 1-N for cols of per-atom array + + int ncompute; // # of Computes accessed by dump + char **id_compute; // their IDs + class Compute **compute; // list of ptrs to the Computes + + int nfix; // # of Fixes used by dump + char **id_fix; // their IDs + class Fix **fix; // list of ptrs to the Fixes + + // private methods + + void init_style() override; + void write_header(bigint) override; + int count() override; + void pack(tagint *) override; + int convert_string(int, double *) override; + void write_data(int, double *) override; + double memory_usage() override; + + int parse_fields(int, char **); + int add_compute(const char *); + int add_fix(const char *); + int modify_param(int, char **) override; + + void header_format_binary(); + void header_unit_style_binary(); + void header_time_binary(); + void header_columns_binary(); + void format_magic_string_binary(); + void format_endian_binary(); + void format_revision_binary(); + + typedef void (DumpGrid::*FnPtrHeader)(bigint); + FnPtrHeader header_choice; // ptr to write header functions + void header_binary(bigint); + void header_binary_triclinic(bigint); + void header_item(bigint); + void header_item_triclinic(bigint); + + typedef void (DumpGrid::*FnPtrWrite)(int, double *); + FnPtrWrite write_choice; // ptr to write data functions + void write_binary(int, double *); + void write_string(int, double *); + void write_lines(int, double *); + + // customize by adding a method prototype + + typedef void (DumpGrid::*FnPtrPack)(int); + FnPtrPack *pack_choice; // ptrs to pack functions + + void pack_compute(int); + void pack_fix(int); +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/fix.h b/src/fix.h index 41a30d2d60..0b15bacff6 100644 --- a/src/fix.h +++ b/src/fix.h @@ -99,6 +99,8 @@ class Fix : protected Pointers { int size_local_cols; // 0 = vector, N = columns in local array int local_freq; // frequency local data is available at + int pergrid_flag; // 0/1 if per-grid data is stored + int extscalar; // 0/1 if global scalar is intensive/extensive int extvector; // 0/1/-1 if global vector is all int/ext/extlist int *extlist; // list of 0/1 int/ext for each vec component diff --git a/src/grid3d.cpp b/src/grid3d.cpp index 8f65cb55b9..f96770d3c9 100644 --- a/src/grid3d.cpp +++ b/src/grid3d.cpp @@ -12,7 +12,7 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#include "grid2d.h" +#include "grid3d.h" #include "comm.h" #include "error.h" @@ -48,7 +48,7 @@ enum{REGULAR,TILED}; communication is done across the periodic boundaries ------------------------------------------------------------------------- */ -Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm, +Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm gcomm, int gnx, int gny, int gnz, int ixlo, int ixhi, int iylo, int iyhi, int izlo, int izhi, int oxlo, int oxhi, int oylo, int oyhi, int ozlo, int ozhi) @@ -87,7 +87,7 @@ Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm, e xyz lohi for flag = 2: 6 neighbor procs ------------------------------------------------------------------------- */ -Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm, int flag, +Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm gcomm, int flag, int gnx, int gny, int gnz, int ixlo, int ixhi, int iylo, int iyhi, int izlo, int izhi, int oxlo, int oxhi, int oylo, int oyhi, int ozlo, int ozhi, @@ -124,14 +124,14 @@ Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm, int flag, oxlo,oxhi,oylo,oyhi,ozlo,ozhi, exlo,exhi,eylo,eyhi,ezlo,ezhi); } else { - error->all(FLERR,"Grid2d does not support tiled layout with neighbor procs"); + error->all(FLERR,"Grid3d does not support tiled layout with neighbor procs"); } } } /* ---------------------------------------------------------------------- */ -Grid2d::~Grid2d() +Grid3d::~Grid3d() { // regular comm data struct @@ -164,7 +164,7 @@ Grid2d::~Grid2d() store constructor args in local variables ------------------------------------------------------------------------- */ -void Grid2d::initialize(MPI_Comm gcomm, +void Grid3d::initialize(MPI_Comm gcomm, int gnx, int gny, int gnz, int ixlo, int ixhi, int iylo, int iyhi, int izlo, int izhi, @@ -229,7 +229,7 @@ void Grid2d::initialize(MPI_Comm gcomm, /* ---------------------------------------------------------------------- */ -void Grid2d::setup(int &nbuf1, int &nbuf2) +void Grid3d::setup(int &nbuf1, int &nbuf2) { if (layout == REGULAR) setup_regular(nbuf1,nbuf2); else setup_tiled(nbuf1,nbuf2); @@ -244,7 +244,7 @@ void Grid2d::setup(int &nbuf1, int &nbuf2) all procs perform same # of swaps in a direction, even if some don't need it ------------------------------------------------------------------------- */ -void Grid2d::setup_regular(int &nbuf1, int &nbuf2) +void Grid3d::setup_regular(int &nbuf1, int &nbuf2) { int nsent,sendfirst,sendlast,recvfirst,recvlast; int sendplanes,recvplanes; @@ -545,7 +545,7 @@ void Grid2d::setup_regular(int &nbuf1, int &nbuf2) no exchanges by dimension, unlike CommTiled forward/reverse comm of particles ------------------------------------------------------------------------- */ -void Grid2d::setup_tiled(int &nbuf1, int &nbuf2) +void Grid3d::setup_tiled(int &nbuf1, int &nbuf2) { int i,m; double xlo,xhi,ylo,yhi,zlo,zhi; @@ -557,7 +557,7 @@ void Grid2d::setup_tiled(int &nbuf1, int &nbuf2) // dim is -1 for proc 0, but never accessed rcbinfo = (RCBinfo *) - memory->smalloc(nprocs*sizeof(RCBinfo),"grid2d:rcbinfo"); + memory->smalloc(nprocs*sizeof(RCBinfo),"grid3d:rcbinfo"); RCBinfo rcbone; rcbone.dim = comm->rcbcutdim; if (rcbone.dim <= 0) rcbone.cut = inxlo; @@ -580,7 +580,7 @@ void Grid2d::setup_tiled(int &nbuf1, int &nbuf2) pbc[0] = pbc[1] = pbc[2] = 0; - memory->create(overlap_procs,nprocs,"grid2d:overlap_procs"); + memory->create(overlap_procs,nprocs,"grid3d:overlap_procs"); noverlap = maxoverlap = 0; overlap = nullptr; @@ -591,9 +591,9 @@ void Grid2d::setup_tiled(int &nbuf1, int &nbuf2) // ncopy = # of overlaps with myself, across a periodic boundary int *proclist; - memory->create(proclist,noverlap,"grid2d:proclist"); + memory->create(proclist,noverlap,"grid3d:proclist"); srequest = (Request *) - memory->smalloc(noverlap*sizeof(Request),"grid2d:srequest"); + memory->smalloc(noverlap*sizeof(Request),"grid3d:srequest"); int nsend_request = 0; ncopy = 0; @@ -612,17 +612,17 @@ void Grid2d::setup_tiled(int &nbuf1, int &nbuf2) auto irregular = new Irregular(lmp); int nrecv_request = irregular->create_data(nsend_request,proclist,1); - auto rrequest = (Request *) memory->smalloc(nrecv_request*sizeof(Request),"grid2d:rrequest"); + auto rrequest = (Request *) memory->smalloc(nrecv_request*sizeof(Request),"grid3d:rrequest"); irregular->exchange_data((char *) srequest,sizeof(Request),(char *) rrequest); irregular->destroy_data(); // compute overlaps between received ghost boxes and my owned box // overlap box used to setup my Send data struct and respond to requests - send = (Send *) memory->smalloc(nrecv_request*sizeof(Send),"grid2d:send"); - sresponse = (Response *) memory->smalloc(nrecv_request*sizeof(Response),"grid2d:sresponse"); + send = (Send *) memory->smalloc(nrecv_request*sizeof(Send),"grid3d:send"); + sresponse = (Response *) memory->smalloc(nrecv_request*sizeof(Response),"grid3d:sresponse"); memory->destroy(proclist); - memory->create(proclist,nrecv_request,"grid2d:proclist"); + memory->create(proclist,nrecv_request,"grid3d:proclist"); for (m = 0; m < nrecv_request; m++) { send[m].proc = rrequest[m].sender; @@ -651,7 +651,7 @@ void Grid2d::setup_tiled(int &nbuf1, int &nbuf2) int nsend_response = nrecv_request; int nrecv_response = irregular->create_data(nsend_response,proclist,1); - auto rresponse = (Response *) memory->smalloc(nrecv_response*sizeof(Response),"grid2d:rresponse"); + auto rresponse = (Response *) memory->smalloc(nrecv_response*sizeof(Response),"grid3d:rresponse"); irregular->exchange_data((char *) sresponse,sizeof(Response),(char *) rresponse); irregular->destroy_data(); delete irregular; @@ -660,7 +660,7 @@ void Grid2d::setup_tiled(int &nbuf1, int &nbuf2) // box used to setup my Recv data struct after unwrapping via PBC // adjacent = 0 if any box of ghost cells does not adjoin my owned cells - recv = (Recv *) memory->smalloc(nrecv_response*sizeof(Recv),"grid2d:recv"); + recv = (Recv *) memory->smalloc(nrecv_response*sizeof(Recv),"grid3d:recv"); adjacent = 1; for (i = 0; i < nrecv_response; i++) { @@ -683,7 +683,7 @@ void Grid2d::setup_tiled(int &nbuf1, int &nbuf2) // create Copy data struct from overlaps with self - copy = (Copy *) memory->smalloc(ncopy*sizeof(Copy),"grid2d:copy"); + copy = (Copy *) memory->smalloc(ncopy*sizeof(Copy),"grid3d:copy"); ncopy = 0; for (m = 0; m < noverlap; m++) { @@ -770,7 +770,7 @@ void Grid2d::setup_tiled(int &nbuf1, int &nbuf2) add all the procs it overlaps with to Overlap list ------------------------------------------------------------------------- */ -void Grid2d::ghost_box_drop(int *box, int *pbc) +void Grid3d::ghost_box_drop(int *box, int *pbc) { int i,m; @@ -855,7 +855,7 @@ void Grid2d::ghost_box_drop(int *box, int *pbc) return Np = # of procs, plist = proc IDs ------------------------------------------------------------------------- */ -void Grid2d::box_drop_grid(int *box, int proclower, int procupper, +void Grid3d::box_drop_grid(int *box, int proclower, int procupper, int &np, int *plist) { // end recursion when partition is a single proc @@ -886,7 +886,7 @@ void Grid2d::box_drop_grid(int *box, int proclower, int procupper, return 1 if yes, 0 if no ------------------------------------------------------------------------- */ -int Grid2d::ghost_adjacent() +int Grid3d::ghost_adjacent() { if (layout == REGULAR) return ghost_adjacent_regular(); return ghost_adjacent_tiled(); @@ -897,7 +897,7 @@ int Grid2d::ghost_adjacent() return 0 if adjacent=0 for any proc, else 1 ------------------------------------------------------------------------- */ -int Grid2d::ghost_adjacent_regular() +int Grid3d::ghost_adjacent_regular() { adjacent = 1; if (ghostxlo > inxhi-inxlo+1) adjacent = 0; @@ -918,7 +918,7 @@ int Grid2d::ghost_adjacent_regular() return 0 if adjacent=0 for any proc, else 1 ------------------------------------------------------------------------- */ -int Grid2d::ghost_adjacent_tiled() +int Grid3d::ghost_adjacent_tiled() { int adjacent_all; MPI_Allreduce(&adjacent,&adjacent_all,1,MPI_INT,MPI_MIN,gridcomm); @@ -929,7 +929,7 @@ int Grid2d::ghost_adjacent_tiled() forward comm of my owned cells to other's ghost cells ------------------------------------------------------------------------- */ -void Grid2d::forward_comm(int caller, void *ptr, int nper, int nbyte, int which, +void Grid3d::forward_comm(int caller, void *ptr, int nper, int nbyte, int which, void *buf1, void *buf2, MPI_Datatype datatype) { if (layout == REGULAR) { @@ -960,7 +960,7 @@ void Grid2d::forward_comm(int caller, void *ptr, int nper, int nbyte, int which, ------------------------------------------------------------------------- */ template < class T > -void Grid2d:: +void Grid3d:: forward_comm_regular(T *ptr, int nper, int /*nbyte*/, int which, void *buf1, void *buf2, MPI_Datatype datatype) { @@ -990,7 +990,7 @@ forward_comm_regular(T *ptr, int nper, int /*nbyte*/, int which, ------------------------------------------------------------------------- */ template < class T > -void Grid2d:: +void Grid3d:: forward_comm_tiled(T *ptr, int nper, int nbyte, int which, void *buf1, void *vbuf2, MPI_Datatype datatype) { @@ -1034,7 +1034,7 @@ forward_comm_tiled(T *ptr, int nper, int nbyte, int which, reverse comm of my ghost cells to sum to owner cells ------------------------------------------------------------------------- */ -void Grid2d::reverse_comm(int caller, void *ptr, int nper, int nbyte, int which, +void Grid3d::reverse_comm(int caller, void *ptr, int nper, int nbyte, int which, void *buf1, void *buf2, MPI_Datatype datatype) { if (layout == REGULAR) { @@ -1065,7 +1065,7 @@ void Grid2d::reverse_comm(int caller, void *ptr, int nper, int nbyte, int which, ------------------------------------------------------------------------- */ template < class T > -void Grid2d:: +void Grid3d:: reverse_comm_regular(T *ptr, int nper, int /*nbyte*/, int which, void *buf1, void *buf2, MPI_Datatype datatype) { @@ -1095,7 +1095,7 @@ reverse_comm_regular(T *ptr, int nper, int /*nbyte*/, int which, ------------------------------------------------------------------------- */ template < class T > -void Grid2d:: +void Grid3d:: reverse_comm_tiled(T *ptr, int nper, int nbyte, int which, void *buf1, void *vbuf2, MPI_Datatype datatype) { @@ -1142,7 +1142,7 @@ reverse_comm_tiled(T *ptr, int nper, int nbyte, int which, caller can decide whether to store chunks, output them, etc ------------------------------------------------------------------------- */ -void Grid2d::gather(int /*caller*/, void *ptr, int nper, int nbyte, +void Grid3d::gather(int /*caller*/, void *ptr, int nper, int nbyte, int which, void *buf, MPI_Datatype datatype) { int me = comm->me; @@ -1158,8 +1158,8 @@ void Grid2d::gather(int /*caller*/, void *ptr, int nper, int nbyte, // pack my data via callback to caller char *mybuf; - if (me == 0) memory->create(mybuf,maxsize*nbyte,"grid2d:mybuf"); - else memory->create(mybuf,mysize*nbyte,"grid2d:mybuf"); + if (me == 0) memory->create(mybuf,maxsize*nbyte,"grid3d:mybuf"); + else memory->create(mybuf,mysize*nbyte,"grid3d:mybuf"); fptr->pack_gather_grid(which,mybuf); // ping each proc for its data @@ -1219,10 +1219,10 @@ void Grid2d::gather(int /*caller*/, void *ptr, int nper, int nbyte, same swap list used by forward and reverse communication ------------------------------------------------------------------------- */ -void Grid2d::grow_swap() +void Grid3d::grow_swap() { maxswap += DELTA; - swap = (Swap *) memory->srealloc(swap,maxswap*sizeof(Swap),"grid2d:swap"); + swap = (Swap *) memory->srealloc(swap,maxswap*sizeof(Swap),"grid3d:swap"); } /* ---------------------------------------------------------------------- @@ -1233,11 +1233,11 @@ void Grid2d::grow_swap() same swap list used by forward and reverse communication ------------------------------------------------------------------------- */ -void Grid2d::grow_overlap() +void Grid3d::grow_overlap() { maxoverlap += DELTA; overlap = (Overlap *) - memory->srealloc(overlap,maxoverlap*sizeof(Overlap),"grid2d:overlap"); + memory->srealloc(overlap,maxoverlap*sizeof(Overlap),"grid3d:overlap"); } /* ---------------------------------------------------------------------- @@ -1246,11 +1246,11 @@ void Grid2d::grow_overlap() (fullxlo:fullxhi,fullylo:fullyhi,fullzlo:fullzhi) ------------------------------------------------------------------------- */ -int Grid2d::indices(int *&list, +int Grid3d::indices(int *&list, int xlo, int xhi, int ylo, int yhi, int zlo, int zhi) { int nmax = (xhi-xlo+1) * (yhi-ylo+1) * (zhi-zlo+1); - memory->create(list,nmax,"grid2d:indices"); + memory->create(list,nmax,"grid3d:indices"); if (nmax == 0) return 0; int nx = (fullxhi-fullxlo+1);