first version of fix ave/grid
This commit is contained in:
Steve Plimpton
@ -291,7 +291,7 @@ void DumpGrid::init_style()
|
||||
}
|
||||
|
||||
} else {
|
||||
if (field2source[0] == COMPUTE) {
|
||||
if (field2source[i] == COMPUTE) {
|
||||
icompute = compute[field2index[i]];
|
||||
grid3d = (Grid3d *) icompute->get_grid_by_index(field2grid[i]);
|
||||
} else {
|
||||
@ -760,10 +760,10 @@ int DumpGrid::parse_fields(int narg, char **arg)
|
||||
if (argi.get_dim() == 0 && ncol)
|
||||
error->all(FLERR,"Dump grid compute {} data {} is not per-grid vector",
|
||||
idcompute,dname);
|
||||
if (argi.get_dim() > 0 && ncol == 0)
|
||||
if (argi.get_dim() && ncol == 0)
|
||||
error->all(FLERR,"Dump grid compute {} data {} is not per-grid array",
|
||||
idcompute,dname);
|
||||
if (argi.get_dim() > 0 && argi.get_index1() > ncol)
|
||||
if (argi.get_dim() && argi.get_index1() > ncol)
|
||||
error->all(FLERR,
|
||||
"Dump grid compute {} array {} is accessed out-of-range",
|
||||
idcompute,dname);
|
||||
@ -798,6 +798,8 @@ int DumpGrid::parse_fields(int narg, char **arg)
|
||||
if (ifix->pergrid_flag == 0)
|
||||
error->all(FLERR,"Dump grid fix {} does not compute per-grid info",
|
||||
idfix);
|
||||
if (update->ntimestep % ifix->pergrid_freq)
|
||||
error->all(FLERR,"Fix for dump grid not computed at compatible time");
|
||||
|
||||
int dim;
|
||||
int igrid = ifix->get_grid_by_name(gname,dim);
|
||||
@ -961,7 +963,7 @@ int DumpGrid::modify_param(int narg, char **arg)
|
||||
double DumpGrid::memory_usage()
|
||||
{
|
||||
double bytes = Dump::memory_usage();
|
||||
//NOTE: restre if use choose
|
||||
//NOTE: restore if use choose
|
||||
//bytes += memory->usage(choose,maxlocal);
|
||||
//bytes += memory->usage(dchoose,maxlocal);
|
||||
//bytes += memory->usage(clist,maxlocal);
|
||||
|
||||
@ -100,6 +100,7 @@ class Fix : protected Pointers {
|
||||
int local_freq; // frequency local data is available at
|
||||
|
||||
int pergrid_flag; // 0/1 if per-grid data is stored
|
||||
int pergrid_freq; // frequency per-grid data is available at
|
||||
|
||||
int extscalar; // 0/1 if global scalar is intensive/extensive
|
||||
int extvector; // 0/1/-1 if global vector is all int/ext/extlist
|
||||
|
||||
765
src/fix_ave_grid.cpp
Normal file
765
src/fix_ave_grid.cpp
Normal file
@ -0,0 +1,765 @@
|
||||
// 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 "fix_ave_grid.h"
|
||||
|
||||
#include "arg_info.h"
|
||||
#include "atom.h"
|
||||
#include "compute.h"
|
||||
#include "domain.h"
|
||||
#include "error.h"
|
||||
#include "grid2d.h"
|
||||
#include "grid3d.h"
|
||||
#include "input.h"
|
||||
#include "memory.h"
|
||||
#include "modify.h"
|
||||
#include "update.h"
|
||||
#include "variable.h"
|
||||
|
||||
#include <cstring>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
enum{SAMPLE,ALL};
|
||||
enum{NOSCALE,ATOM};
|
||||
enum{ONE,RUNNING,WINDOW};
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixAveGrid::FixAveGrid(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg),
|
||||
which(nullptr), argindex(nullptr), ids(nullptr),
|
||||
value2index(nullptr), value2grid(nullptr), value2data(nullptr),
|
||||
grid2d(nullptr), grid3d(nullptr),
|
||||
vec2d(nullptr), array2d(nullptr), vec3d(nullptr), array3d(nullptr)
|
||||
{
|
||||
if (narg < 10) error->all(FLERR,"Illegal fix ave/grid command");
|
||||
|
||||
pergrid_flag = 1;
|
||||
nevery = utils::inumeric(FLERR,arg[3],false,lmp);
|
||||
nrepeat = utils::inumeric(FLERR,arg[4],false,lmp);
|
||||
pergrid_freq = utils::inumeric(FLERR,arg[5],false,lmp);
|
||||
time_depend = 1;
|
||||
|
||||
// NOTE: test for Dxyz as well
|
||||
|
||||
nx = utils::inumeric(FLERR,arg[6],false,lmp);
|
||||
ny = utils::inumeric(FLERR,arg[7],false,lmp);
|
||||
nz = utils::inumeric(FLERR,arg[8],false,lmp);
|
||||
|
||||
// expand args if any have wildcard character "*"
|
||||
// this can reset nvalues
|
||||
|
||||
int expand = 0;
|
||||
char **earg;
|
||||
int nargnew = utils::expand_args(FLERR,nvalues,&arg[9],1,earg,lmp);
|
||||
|
||||
if (earg != &arg[9]) expand = 1;
|
||||
arg = earg;
|
||||
|
||||
// parse values
|
||||
|
||||
which = new int[nvalues];
|
||||
argindex = new int[nvalues];
|
||||
ids = new char*[nvalues];
|
||||
value2index = new int[nvalues];
|
||||
value2grid = new int[nvalues];
|
||||
value2data = new int[nvalues];
|
||||
|
||||
modeatom = modegrid = 0;
|
||||
|
||||
int iarg = 0;
|
||||
while (iarg < nargnew) {
|
||||
ids[nvalues] = nullptr;
|
||||
|
||||
if (strcmp(arg[iarg],"vx") == 0) {
|
||||
which[nvalues] = ArgInfo::V;
|
||||
argindex[nvalues++] = 0;
|
||||
modeatom = 1;
|
||||
} else if (strcmp(arg[iarg],"vy") == 0) {
|
||||
which[nvalues] = ArgInfo::V;
|
||||
argindex[nvalues++] = 1;
|
||||
modeatom = 1;
|
||||
} else if (strcmp(arg[iarg],"vz") == 0) {
|
||||
which[nvalues] = ArgInfo::V;
|
||||
argindex[nvalues++] = 2;
|
||||
modeatom = 1;
|
||||
|
||||
} else if (strcmp(arg[iarg],"fx") == 0) {
|
||||
which[nvalues] = ArgInfo::F;
|
||||
argindex[nvalues++] = 0;
|
||||
modeatom = 1;
|
||||
} else if (strcmp(arg[iarg],"fy") == 0) {
|
||||
which[nvalues] = ArgInfo::F;
|
||||
argindex[nvalues++] = 1;
|
||||
modeatom = 1;
|
||||
} else if (strcmp(arg[iarg],"fz") == 0) {
|
||||
which[nvalues] = ArgInfo::F;
|
||||
argindex[nvalues++] = 2;
|
||||
modeatom = 1;
|
||||
|
||||
} else if (strcmp(arg[iarg],"density/number") == 0) {
|
||||
which[nvalues] = ArgInfo::DENSITY_NUMBER;
|
||||
argindex[nvalues++] = 0;
|
||||
modeatom = 1;
|
||||
} else if (strcmp(arg[iarg],"density/mass") == 0) {
|
||||
which[nvalues] = ArgInfo::DENSITY_MASS;
|
||||
argindex[nvalues++] = 0;
|
||||
modeatom = 1;
|
||||
} else if (strcmp(arg[iarg],"mass") == 0) {
|
||||
which[nvalues] = ArgInfo::MASS;
|
||||
argindex[nvalues++] = 0;
|
||||
modeatom = 1;
|
||||
} else if (strcmp(arg[iarg],"temp") == 0) {
|
||||
which[nvalues] = ArgInfo::TEMPERATURE;
|
||||
argindex[nvalues++] = 0;
|
||||
modeatom = 1;
|
||||
|
||||
} else {
|
||||
ArgInfo argi(arg[iarg]);
|
||||
|
||||
if (argi.get_type() == ArgInfo::NONE) break;
|
||||
if ((argi.get_type() == ArgInfo::UNKNOWN) || (argi.get_dim() > 1))
|
||||
error->all(FLERR,"Invalid fix ave/grid command");
|
||||
|
||||
which[nvalues] = argi.get_type();
|
||||
argindex[nvalues] = argi.get_index1();
|
||||
ids[nvalues] = argi.copy_name();
|
||||
|
||||
if (strchr(ids[nvalues],':')) modegrid = 1;
|
||||
else modeatom = 1;
|
||||
|
||||
if (modegrid && which[nvalues] == ArgInfo::VARIABLE)
|
||||
error->all(FLERR,"Fix ave/grid cannot use variable for grid info");
|
||||
|
||||
nvalues++;
|
||||
}
|
||||
|
||||
iarg++;
|
||||
}
|
||||
|
||||
if (nvalues == 0) error->all(FLERR,"No values in fix ave/grid command");
|
||||
if (modeatom && modegrid)
|
||||
error->all(FLERR,"Fix ave/grid cannot operate on per-atom and "
|
||||
"per-grid values");
|
||||
|
||||
// optional args
|
||||
|
||||
normflag = ALL;
|
||||
scaleflag = ATOM;
|
||||
ave = ONE;
|
||||
nwindow = 0;
|
||||
|
||||
while (iarg < nargnew) {
|
||||
if (strcmp(arg[iarg],"norm") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/grid command");
|
||||
if (strcmp(arg[iarg+1],"all") == 0) {
|
||||
normflag = ALL;
|
||||
scaleflag = ATOM;
|
||||
} else if (strcmp(arg[iarg+1],"sample") == 0) {
|
||||
normflag = SAMPLE;
|
||||
scaleflag = ATOM;
|
||||
} else if (strcmp(arg[iarg+1],"none") == 0) {
|
||||
normflag = SAMPLE;
|
||||
scaleflag = NOSCALE;
|
||||
} else error->all(FLERR,"Illegal fix ave/grid command");
|
||||
iarg += 2;
|
||||
} else if (strcmp(arg[iarg],"ave") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/grid command");
|
||||
if (strcmp(arg[iarg+1],"one") == 0) ave = ONE;
|
||||
else if (strcmp(arg[iarg+1],"running") == 0) ave = RUNNING;
|
||||
else if (strcmp(arg[iarg+1],"window") == 0) ave = WINDOW;
|
||||
else error->all(FLERR,"Illegal fix ave/grid command");
|
||||
if (ave == WINDOW) {
|
||||
if (iarg+3 > narg) error->all(FLERR,"Illegal fix ave/grid command");
|
||||
nwindow = utils::inumeric(FLERR,arg[iarg+2],false,lmp);
|
||||
if (nwindow <= 0) error->all(FLERR,"Illegal fix ave/grid command");
|
||||
}
|
||||
iarg += 2;
|
||||
if (ave == WINDOW) iarg++;
|
||||
}
|
||||
}
|
||||
|
||||
// if wildcard expansion occurred, free earg memory from exapnd_args()
|
||||
|
||||
if (expand) {
|
||||
for (int i = 0; i < nvalues; i++) delete [] earg[i];
|
||||
memory->sfree(earg);
|
||||
}
|
||||
|
||||
// setup and error check
|
||||
// for fix inputs, check that fix frequency is acceptable
|
||||
|
||||
dimension = domain->dimension;
|
||||
|
||||
if (nevery <= 0 || nrepeat <= 0 || pergrid_freq <= 0)
|
||||
error->all(FLERR,"Illegal fix ave/grid command");
|
||||
if (pergrid_freq % nevery || nrepeat*nevery > pergrid_freq)
|
||||
error->all(FLERR,"Illegal fix ave/grid command");
|
||||
|
||||
if (nx < 1 || ny < 1 || nz < 1)
|
||||
error->all(FLERR,"Invalid fix ave/grid grid size");
|
||||
if (dimension == 2 && nz != 1)
|
||||
error->all(FLERR,"Fix ave/grid grid Nz must be 1 for 2d simulation");
|
||||
|
||||
// error checks for ATOM mode
|
||||
|
||||
if (modeatom) {
|
||||
for (int i = 0; i < nvalues; i++) {
|
||||
if (which[i] == ArgInfo::COMPUTE) {
|
||||
int icompute = modify->find_compute(ids[i]);
|
||||
if (icompute < 0)
|
||||
error->all(FLERR,"Compute ID for fix ave/grid does not exist");
|
||||
if (modify->compute[icompute]->peratom_flag == 0)
|
||||
error->all(FLERR,
|
||||
"Fix ave/atom compute does not calculate per-atom values");
|
||||
if (argindex[i] == 0 &&
|
||||
modify->compute[icompute]->size_peratom_cols != 0)
|
||||
error->all(FLERR,"Fix ave/atom compute does not "
|
||||
"calculate a per-atom vector");
|
||||
if (argindex[i] && modify->compute[icompute]->size_peratom_cols == 0)
|
||||
error->all(FLERR,"Fix ave/atom compute does not "
|
||||
"calculate a per-atom array");
|
||||
if (argindex[i] &&
|
||||
argindex[i] > modify->compute[icompute]->size_peratom_cols)
|
||||
error->all(FLERR,"Fix ave/atom compute array is accessed out-of-range");
|
||||
|
||||
} else if (which[i] == ArgInfo::FIX) {
|
||||
int ifix = modify->find_fix(ids[i]);
|
||||
if (ifix < 0)
|
||||
error->all(FLERR,"Fix ID for fix ave/atom does not exist");
|
||||
if (modify->fix[ifix]->peratom_flag == 0)
|
||||
error->all(FLERR,"Fix ave/atom fix does not calculate per-atom values");
|
||||
if (argindex[i] == 0 && modify->fix[ifix]->size_peratom_cols != 0)
|
||||
error->all(FLERR,
|
||||
"Fix ave/atom fix does not calculate a per-atom vector");
|
||||
if (argindex[i] && modify->fix[ifix]->size_peratom_cols == 0)
|
||||
error->all(FLERR,
|
||||
"Fix ave/atom fix does not calculate a per-atom array");
|
||||
if (argindex[i] && argindex[i] > modify->fix[ifix]->size_peratom_cols)
|
||||
error->all(FLERR,"Fix ave/atom fix array is accessed out-of-range");
|
||||
if (nevery % modify->fix[ifix]->peratom_freq)
|
||||
error->all(FLERR,
|
||||
"Fix for fix ave/atom not computed at compatible time");
|
||||
|
||||
} else if (which[i] == ArgInfo::VARIABLE) {
|
||||
int ivariable = input->variable->find(ids[i]);
|
||||
if (ivariable < 0)
|
||||
error->all(FLERR,"Variable name for fix ave/atom does not exist");
|
||||
if (input->variable->atomstyle(ivariable) == 0)
|
||||
error->all(FLERR,"Fix ave/atom variable is not atom-style variable");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// setup and error checks for GRID mode
|
||||
|
||||
if (modegrid) {
|
||||
for (int i = 0; i < nvalues; i++) {
|
||||
if (which[i] == ArgInfo::COMPUTE) {
|
||||
|
||||
char *idcompute,*gname,*dname;
|
||||
utils::grid_parse(FLERR,ids[i],idcompute,gname,dname,error);
|
||||
|
||||
Compute *icompute = modify->get_compute_by_id(idcompute);
|
||||
if (!icompute)
|
||||
error->all(FLERR,"Could not find fix ave/grid compute ID: {}",
|
||||
idcompute);
|
||||
if (icompute->pergrid_flag == 0)
|
||||
error->all(FLERR,
|
||||
"Fix ave/grid compute {} does not compute per-grid info",
|
||||
idcompute);
|
||||
|
||||
int dim;
|
||||
int igrid = icompute->get_grid_by_name(gname,dim);
|
||||
if (igrid < 0)
|
||||
error->all(FLERR,
|
||||
"Fix ave/grid compute {} does not recognize grid name {}",
|
||||
idcompute,gname);
|
||||
|
||||
int ncol;
|
||||
int idata = icompute->get_griddata_by_name(igrid,dname,ncol);
|
||||
if (idata < 0)
|
||||
error->all(FLERR,
|
||||
"Fix ave/grid compute {} does not recognize data name {}",
|
||||
idcompute,dname);
|
||||
|
||||
if (argindex[i] == 0 && ncol)
|
||||
error->all(FLERR,
|
||||
"Fix ave/grid compute {} data {} is not per-grid vector",
|
||||
idcompute,dname);
|
||||
if (argindex[i] && ncol == 0)
|
||||
error->all(FLERR,
|
||||
"Fix ave/grid compute {} data {} is not per-grid array",
|
||||
idcompute,dname);
|
||||
if (argindex[i] && argindex[i] > ncol)
|
||||
error->all(FLERR,
|
||||
"Fix ave/grid compute {} array {} is accessed out-of-range",
|
||||
idcompute,dname);
|
||||
|
||||
value2grid[iarg] = igrid;
|
||||
value2data[iarg] = idata;
|
||||
|
||||
delete [] idcompute;
|
||||
delete [] gname;
|
||||
delete [] dname;
|
||||
|
||||
} else if (which[i] == ArgInfo::FIX) {
|
||||
|
||||
char *idfix,*gname,*dname;
|
||||
utils::grid_parse(FLERR,ids[i],idfix,gname,dname,error);
|
||||
|
||||
Fix *ifix = modify->get_fix_by_id(idfix);
|
||||
if (!ifix) error->all(FLERR,"Could not find fix ave/grid fix ID: {}",
|
||||
idfix);
|
||||
if (ifix->pergrid_flag == 0)
|
||||
error->all(FLERR,"Fix ave/grid fix {} does not compute per-grid info",
|
||||
idfix);
|
||||
if (nevery % ifix->pergrid_freq)
|
||||
error->all(FLERR,
|
||||
"Fix for fix grid/atom not computed at compatible time");
|
||||
|
||||
int dim;
|
||||
int igrid = ifix->get_grid_by_name(gname,dim);
|
||||
if (igrid < 0)
|
||||
error->all(FLERR,
|
||||
"Fix ave/grid compute {} does not recognize grid name {}",
|
||||
idfix,gname);
|
||||
|
||||
int ncol;
|
||||
int idata = ifix->get_griddata_by_name(igrid,dname,ncol);
|
||||
if (idata < 0)
|
||||
error->all(FLERR,
|
||||
"Fix ave/grid compute {} does not recognize data name {}",
|
||||
idfix,dname);
|
||||
|
||||
if (argindex[i] == 0 && ncol)
|
||||
error->all(FLERR,
|
||||
"Fix ave/grid compute {} data {} is not per-grid vector",
|
||||
idfix,dname);
|
||||
if (argindex[i] && ncol == 0)
|
||||
error->all(FLERR,
|
||||
"Fix ave/grid compute {} data {} is not per-grid array",
|
||||
idfix,dname);
|
||||
if (argindex[i] && argindex[i] > ncol)
|
||||
error->all(FLERR,
|
||||
"Fix ave/grid compute {} array {} is accessed out-of-range",
|
||||
idfix,dname);
|
||||
|
||||
value2grid[iarg] = igrid;
|
||||
value2data[iarg] = idata;
|
||||
|
||||
delete [] idfix;
|
||||
delete [] gname;
|
||||
delete [] dname;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// instantiate the Grid class and allocate per-grid memory
|
||||
// NOTE: need to extend ghost grid for ATOM mode ?
|
||||
|
||||
if (dimension == 2) {
|
||||
grid2d = new Grid2d(lmp, world, nx, ny, 0, 0.0, 0.0,
|
||||
nxlo_in, nxhi_in, nylo_in, nyhi_in,
|
||||
nxlo_out, nxhi_out, nylo_out, nyhi_out);
|
||||
if (nvalues == 1)
|
||||
memory->create2d_offset(vec2d, nylo_out, nyhi_out, nxlo_out, nxhi_out,
|
||||
"fix_ave/grid:vec2d");
|
||||
else
|
||||
memory->create3d_offset_last(array2d, nylo_out, nyhi_out, nxlo_out,
|
||||
nxhi_out, nvalues, "fix_ave/grid:array2d");
|
||||
ngridout = (nxhi_out - nxlo_out + 1) * (nyhi_out - nylo_out + 1);
|
||||
|
||||
} else {
|
||||
grid3d = new Grid3d(lmp, world, nx, ny, nz, 0, 0.0, 0.0,
|
||||
nxlo_in, nxhi_in, nylo_in, nyhi_in, nzlo_in, nzhi_in,
|
||||
nxlo_out, nxhi_out, nylo_out, nyhi_out,
|
||||
nzlo_out, nzhi_out);
|
||||
if (nvalues == 1)
|
||||
memory->create3d_offset(vec3d, nzlo_out, nzhi_out, nylo_out,
|
||||
nyhi_out, nxlo_out,
|
||||
nxhi_out, "fix_ave/grid:vec3d");
|
||||
else
|
||||
memory->create4d_offset_last(array3d, nzlo_out, nzhi_out, nylo_out,
|
||||
nyhi_out, nxlo_out,
|
||||
nxhi_out, nvalues, "fix_ave/grid:array3d");
|
||||
ngridout = (nxhi_out - nxlo_out + 1) * (nyhi_out - nylo_out + 1) *
|
||||
(nzhi_out - nzlo_out + 1);
|
||||
}
|
||||
|
||||
// zero the array since dump may access it on timestep 0
|
||||
// zero the array since a variable may access it before first run
|
||||
|
||||
/*
|
||||
int nlocal = atom->nlocal;
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
for (int m = 0; m < nvalues; m++)
|
||||
array[i][m] = 0.0;
|
||||
*/
|
||||
|
||||
// nvalid = next step on which end_of_step does something
|
||||
// add nvalid to all computes that store invocation times
|
||||
// since don't know a priori which are invoked by this fix
|
||||
// once in end_of_step() can set timestep for ones actually invoked
|
||||
|
||||
irepeat = 0;
|
||||
nvalid_last = -1;
|
||||
nvalid = nextvalid();
|
||||
modify->addstep_compute_all(nvalid);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixAveGrid::~FixAveGrid()
|
||||
{
|
||||
delete [] which;
|
||||
delete [] argindex;
|
||||
for (int m = 0; m < nvalues; m++) delete [] ids[m];
|
||||
delete [] ids;
|
||||
delete [] value2index;
|
||||
delete [] value2grid;
|
||||
delete [] value2data;
|
||||
|
||||
delete grid2d;
|
||||
delete grid3d;
|
||||
memory->destroy2d_offset(vec2d,nylo_out,nxlo_out);
|
||||
memory->destroy2d_offset(array2d,nylo_out,nxlo_out);
|
||||
memory->destroy3d_offset(vec3d,nzlo_out,nylo_out,nxlo_out);
|
||||
memory->destroy4d_offset_last(array3d,nzlo_out,nylo_out,nxlo_out);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixAveGrid::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= END_OF_STEP;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixAveGrid::init()
|
||||
{
|
||||
// set indices and check validity of all computes,fixes,variables
|
||||
|
||||
for (int m = 0; m < nvalues; m++) {
|
||||
if (which[m] == ArgInfo::COMPUTE) {
|
||||
int icompute = modify->find_compute(ids[m]);
|
||||
if (icompute < 0)
|
||||
error->all(FLERR,"Compute ID for fix ave/atom does not exist");
|
||||
value2index[m] = icompute;
|
||||
|
||||
} else if (which[m] == ArgInfo::FIX) {
|
||||
int ifix = modify->find_fix(ids[m]);
|
||||
if (ifix < 0)
|
||||
error->all(FLERR,"Fix ID for fix ave/atom does not exist");
|
||||
value2index[m] = ifix;
|
||||
|
||||
} else if (which[m] == ArgInfo::VARIABLE) {
|
||||
int ivariable = input->variable->find(ids[m]);
|
||||
if (ivariable < 0)
|
||||
error->all(FLERR,"Variable name for fix ave/atom does not exist");
|
||||
value2index[m] = ivariable;
|
||||
|
||||
} else value2index[m] = -1;
|
||||
}
|
||||
|
||||
// check that grid sizes for all fields match grid size for this fix
|
||||
|
||||
if (modegrid) {
|
||||
Compute *icompute;
|
||||
Fix *ifix;
|
||||
Grid2d *grid2d;
|
||||
Grid3d *grid3d;
|
||||
|
||||
int nxtmp,nytmp,nztmp;
|
||||
|
||||
for (int m = 0; m < nvalues; m++) {
|
||||
if (dimension == 2) {
|
||||
if (which[m] == ArgInfo::COMPUTE) {
|
||||
icompute = modify->compute[value2index[m]];
|
||||
grid2d = (Grid2d *) icompute->get_grid_by_index(value2grid[m]);
|
||||
} else {
|
||||
ifix = modify->fix[value2index[m]];
|
||||
grid2d = (Grid2d *) ifix->get_grid_by_index(value2grid[m]);
|
||||
}
|
||||
grid2d->get_size(nxtmp,nytmp);
|
||||
if (nxtmp != nx || nytmp != ny)
|
||||
error->all(FLERR,"Fix ave/grid value grid sizes do not match");
|
||||
|
||||
} else {
|
||||
if (which[m] == ArgInfo::COMPUTE) {
|
||||
icompute = modify->compute[value2index[m]];
|
||||
grid3d = (Grid3d *) icompute->get_grid_by_index(value2grid[m]);
|
||||
} else {
|
||||
ifix = modify->fix[value2index[m]];
|
||||
grid3d = (Grid3d *) ifix->get_grid_by_index(value2grid[m]);
|
||||
}
|
||||
grid3d->get_size(nxtmp,nytmp,nztmp);
|
||||
if (nxtmp != nx || nytmp != ny || nztmp != nz)
|
||||
error->all(FLERR,"Fix ave/grid value grid sizes do not match");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// need to reset nvalid if nvalid < ntimestep b/c minimize was performed
|
||||
|
||||
if (nvalid < update->ntimestep) {
|
||||
irepeat = 0;
|
||||
nvalid = nextvalid();
|
||||
modify->addstep_compute_all(nvalid);
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
only does something if nvalid = current timestep
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixAveGrid::setup(int /*vflag*/)
|
||||
{
|
||||
end_of_step();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixAveGrid::end_of_step()
|
||||
{
|
||||
int i,j,m,n;
|
||||
|
||||
// skip if not step which requires doing something
|
||||
|
||||
bigint ntimestep = update->ntimestep;
|
||||
if (ntimestep != nvalid) return;
|
||||
nvalid_last = nvalid;
|
||||
|
||||
// zero if first step
|
||||
|
||||
/*
|
||||
if (irepeat == 0)
|
||||
for (i = 0; i < nlocal; i++)
|
||||
for (m = 0; m < nvalues; m++)
|
||||
array[i][m] = 0.0;
|
||||
*/
|
||||
|
||||
// ATOM mode
|
||||
// accumulate per-atom attributes,computes,fixes,variables to local grid
|
||||
// compute/fix/variable may invoke computes so wrap with clear/add
|
||||
|
||||
/*
|
||||
if (modeatom) {
|
||||
modify->clearstep_compute();
|
||||
|
||||
int *mask = atom->mask;
|
||||
|
||||
for (m = 0; m < nvalues; m++) {
|
||||
n = value2index[m];
|
||||
j = argindex[m];
|
||||
|
||||
if (which[m] == ArgInfo::X) {
|
||||
double **x = atom->x;
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) array[i][m] += x[i][j];
|
||||
|
||||
} else if (which[m] == ArgInfo::V) {
|
||||
double **v = atom->v;
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) array[i][m] += v[i][j];
|
||||
|
||||
} else if (which[m] == ArgInfo::F) {
|
||||
double **f = atom->f;
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) array[i][m] += f[i][j];
|
||||
|
||||
// invoke compute if not previously invoked
|
||||
|
||||
} else if (which[m] == ArgInfo::COMPUTE) {
|
||||
Compute *compute = modify->compute[n];
|
||||
if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) {
|
||||
compute->compute_peratom();
|
||||
compute->invoked_flag |= Compute::INVOKED_PERATOM;
|
||||
}
|
||||
|
||||
if (j == 0) {
|
||||
double *compute_vector = compute->vector_atom;
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) array[i][m] += compute_vector[i];
|
||||
} else {
|
||||
int jm1 = j - 1;
|
||||
double **compute_array = compute->array_atom;
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) array[i][m] += compute_array[i][jm1];
|
||||
}
|
||||
|
||||
// access fix fields, guaranteed to be ready
|
||||
|
||||
} else if (which[m] == ArgInfo::FIX) {
|
||||
if (j == 0) {
|
||||
double *fix_vector = modify->fix[n]->vector_atom;
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) array[i][m] += fix_vector[i];
|
||||
} else {
|
||||
int jm1 = j - 1;
|
||||
double **fix_array = modify->fix[n]->array_atom;
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) array[i][m] += fix_array[i][jm1];
|
||||
}
|
||||
|
||||
// evaluate atom-style variable
|
||||
// final argument = 1 sums result to array
|
||||
|
||||
} else if (which[m] == ArgInfo::VARIABLE) {
|
||||
if (array) input->variable->compute_atom(n,igroup,&array[0][m],nvalues,1);
|
||||
else input->variable->compute_atom(n,igroup,nullptr,nvalues,1);
|
||||
}
|
||||
}
|
||||
}
|
||||
*/
|
||||
|
||||
// GRID mode
|
||||
// accumulate results of computes & fixes to local grid
|
||||
|
||||
if (modegrid) {
|
||||
}
|
||||
|
||||
// done if irepeat < nrepeat
|
||||
// else reset irepeat and nvalid
|
||||
|
||||
irepeat++;
|
||||
if (irepeat < nrepeat) {
|
||||
nvalid += nevery;
|
||||
modify->addstep_compute(nvalid);
|
||||
return;
|
||||
}
|
||||
|
||||
irepeat = 0;
|
||||
nvalid = ntimestep+peratom_freq - ((bigint)nrepeat-1)*nevery;
|
||||
modify->addstep_compute(nvalid);
|
||||
|
||||
// if (array == nullptr) return;
|
||||
|
||||
|
||||
//NOTE: need to do comm for atom mode ?
|
||||
|
||||
|
||||
// average the final result for the Nfreq timestep
|
||||
|
||||
/*
|
||||
double repeat = nrepeat;
|
||||
for (i = 0; i < nlocal; i++)
|
||||
for (m = 0; m < nvalues; m++)
|
||||
array[i][m] /= repeat;
|
||||
*/
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
return index of grid associated with name
|
||||
this class can store M named grids, indexed 0 to M-1
|
||||
also set dim for 2d vs 3d grid
|
||||
return -1 if grid name not found
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
int FixAveGrid::get_grid_by_name(char *name, int &dim)
|
||||
{
|
||||
if (strcmp(name,"grid") == 0) {
|
||||
dim = dimension;
|
||||
return 0;
|
||||
}
|
||||
|
||||
return -1;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
return ptr to Grid data struct for grid with index
|
||||
this class can store M named grids, indexed 0 to M-1
|
||||
return nullptr if index is invalid
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void *FixAveGrid::get_grid_by_index(int index)
|
||||
{
|
||||
if (index == 0) {
|
||||
if (dimension == 2) return grid2d;
|
||||
else return grid3d;
|
||||
}
|
||||
return nullptr;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
return index of data associated with name in grid with index igrid
|
||||
this class can store M named grids, indexed 0 to M-1
|
||||
each grid can store G named data sets, indexed 0 to G-1
|
||||
a data set name can be associated with multiple grids
|
||||
set ncol for data set, 0 = vector, 1-N for array with N columns
|
||||
vector = single value per grid pt, array = N values per grid pt
|
||||
return -1 if data name not found
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
int FixAveGrid::get_griddata_by_name(int igrid, char *name, int &ncol)
|
||||
{
|
||||
if (igrid == 0 && strcmp(name,"data") == 0) {
|
||||
if (nvalues == 1) ncol = 0;
|
||||
else ncol = nvalues;
|
||||
return 0;
|
||||
}
|
||||
|
||||
return -1;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
return ptr to multidim data array associated with index
|
||||
this class can store G named data sets, indexed 0 to M-1
|
||||
return nullptr if index is invalid
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void *FixAveGrid::get_griddata_by_index(int index)
|
||||
{
|
||||
if (index == 0) {
|
||||
if (dimension == 2) {
|
||||
if (nvalues == 1) return vec2d;
|
||||
else return array2d;
|
||||
} else {
|
||||
if (nvalues == 1) return vec3d;
|
||||
else return array3d;
|
||||
}
|
||||
}
|
||||
return nullptr;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
memory usage of local atom-based array
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double FixAveGrid::memory_usage()
|
||||
{
|
||||
double bytes = (double) ngridout * nvalues * sizeof(double);
|
||||
return bytes;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
calculate nvalid = next step on which end_of_step does something
|
||||
can be this timestep if multiple of nfreq and nrepeat = 1
|
||||
else backup from next multiple of nfreq
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
bigint FixAveGrid::nextvalid()
|
||||
{
|
||||
bigint nvalid = (update->ntimestep/peratom_freq)*peratom_freq + peratom_freq;
|
||||
if (nvalid-peratom_freq == update->ntimestep && nrepeat == 1)
|
||||
nvalid = update->ntimestep;
|
||||
else
|
||||
nvalid -= ((bigint)nrepeat-1)*nevery;
|
||||
if (nvalid < update->ntimestep) nvalid += peratom_freq;
|
||||
return nvalid;
|
||||
}
|
||||
70
src/fix_ave_grid.h
Normal file
70
src/fix_ave_grid.h
Normal file
@ -0,0 +1,70 @@
|
||||
/* -*- 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 FIX_CLASS
|
||||
// clang-format off
|
||||
FixStyle(ave/grid,FixAveGrid);
|
||||
// clang-format on
|
||||
#else
|
||||
|
||||
#ifndef LMP_FIX_AVE_GRID_H
|
||||
#define LMP_FIX_AVE_GRID_H
|
||||
|
||||
#include "fix.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixAveGrid : public Fix {
|
||||
public:
|
||||
FixAveGrid(class LAMMPS *, int, char **);
|
||||
~FixAveGrid() override;
|
||||
int setmask() override;
|
||||
void init() override;
|
||||
void setup(int) override;
|
||||
void end_of_step() override;
|
||||
int get_grid_by_name(char *, int &) override;
|
||||
void *get_grid_by_index(int) override;
|
||||
int get_griddata_by_name(int, char *, int &) override;
|
||||
void *get_griddata_by_index(int) override;
|
||||
double memory_usage() override;
|
||||
|
||||
private:
|
||||
int nx,ny,nz;
|
||||
int nvalues;
|
||||
int nrepeat, irepeat;
|
||||
bigint nvalid, nvalid_last;
|
||||
int modeatom,modegrid;
|
||||
int normflag,scaleflag,ave,nwindow;
|
||||
int dimension;
|
||||
|
||||
int *which, *argindex;
|
||||
char **ids;
|
||||
int *value2index, *value2grid, *value2data;
|
||||
|
||||
class Grid2d *grid2d;
|
||||
class Grid3d *grid3d;
|
||||
|
||||
int nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in;
|
||||
int nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out;
|
||||
int ngridout;
|
||||
|
||||
double **vec2d,***vec3d;
|
||||
double ***array2d,****array3d;
|
||||
|
||||
bigint nextvalid();
|
||||
};
|
||||
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
Reference in New Issue
Block a user