Prototyping example and accessory commands

This commit is contained in:
jtclemm
2023-04-25 13:38:37 -06:00
parent 7cfe45c00b
commit 7fc916a1d4
5 changed files with 261 additions and 6 deletions

View File

@ -2,6 +2,6 @@ RHEO or Reproducing Hydrodynamics and Elastic Objects is a package to model mult
systems. The authors include Joel Clemmer (Sandia), Thomas O'Connor (Carnegie Mellon), and
Eric Palermo (Carnegie Mellon).
This package requires the GNU scientific library (GSL) version 2.7 or later. To build this
package, one must first separately install GSL in a location that can be found by your
environment.
This package requires the GNU scientific library (GSL). We recommend version 2.7 or later. To
build this package, one must first separately install GSL in a location that can be found by
your environment.

View File

@ -0,0 +1,162 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors:
Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU)
----------------------------------------------------------------------- */
#include "compute_rheo_property_atom.h"
#include "atom.h"
#include "atom_vec.h"
#include "error.h"
#include "memory.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
index(nullptr), colindex(nullptr), pack_choice(nullptr)
{
if (narg < 4) utils::missing_cmd_args(FLERR, "compute property/atom", error);
peratom_flag = 1;
nvalues = narg - 3;
if (nvalues == 1) size_peratom_cols = 0;
else size_peratom_cols = nvalues;
// parse input values
// customize a new keyword by adding to if statement
pack_choice = new FnPtrPack[nvalues];
int i;
for (int iarg = 3; iarg < narg; iarg++) {
i = iarg-3;
if (strcmp(arg[iarg],"id") == 0) {
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_id;
} else if (strcmp(arg[iarg],"mol") == 0) {
if (!atom->molecule_flag)
error->all(FLERR,"Compute property/atom {} is not available", arg[iarg]);
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_molecule;
} else if (strcmp(arg[iarg],"proc") == 0) {
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_proc;
} else if (strcmp(arg[iarg],"type") == 0) {
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_type;
} else if (strcmp(arg[iarg],"mass") == 0) {
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_mass;
} else {
error->all(FLERR,"Invalid keyword {} for compute rheo/property/atom command ", arg[iarg]);
}
}
nmax = 0;
}
/* ---------------------------------------------------------------------- */
ComputeRHEOPropertyAtom::~ComputeRHEOPropertyAtom()
{
delete[] pack_choice;
memory->destroy(vector_atom);
memory->destroy(array_atom);
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::compute_peratom()
{
invoked_peratom = update->ntimestep;
// grow vector or array if necessary
if (atom->nmax > nmax) {
nmax = atom->nmax;
if (nvalues == 1) {
memory->destroy(vector_atom);
memory->create(vector_atom,nmax,"rheo/property/atom:vector");
} else {
memory->destroy(array_atom);
memory->create(array_atom,nmax,nvalues,"rheo/property/atom:array");
}
}
// fill vector or array with per-atom values
if (nvalues == 1) {
buf = vector_atom;
(this->*pack_choice[0])(0);
} else {
if (nmax) buf = &array_atom[0][0];
else buf = nullptr;
for (int n = 0; n < nvalues; n++)
(this->*pack_choice[n])(n);
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeRHEOPropertyAtom::memory_usage()
{
double bytes = (double)nmax * nvalues * sizeof(double);
return bytes;
}
/* ----------------------------------------------------------------------
one method for every keyword compute rheo/property/atom can output
the atom property is packed into buf starting at n with stride nvalues
customize a new keyword by adding a method
------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_rho(int n)
{
double *rho = atom->rho;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = rho[i];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_status(int n)
{
int *status = atom->status;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = status[i];
else buf[n] = 0.0;
n += nvalues;
}
}

View File

@ -0,0 +1,68 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
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 COMPUTE_CLASS
// clang-format off
ComputeStyle(rheo/property/atom,ComputeRHEOPropertyAtom);
// clang-format on
#else
#ifndef LMP_COMPUTE_RHEO_PROPERTY_ATOM_H
#define LMP_COMPUTE_RHEO_PROPERTY_ATOM_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeRHEOPropertyAtom : public Compute {
public:
ComputeRHEOPropertyAtom(class LAMMPS *, int, char **);
~ComputeRHEOPropertyAtom() override;
void compute_peratom() override;
double memory_usage() override;
private:
int nvalues;
int nmax;
double *buf;
typedef void (ComputeRHEOPropertyAtom::*FnPtrPack)(int);
FnPtrPack *pack_choice; // ptrs to pack functions
void pack_rho(int);
void pack_drho(int);
void pack_temperature(int);
void pack_heatflow(int);
void pack_status(int);
void pack_phase(int);
void pack_surface(int);
void pack_r_surface(int);
void pack_divr_surface(int);
void pack_nx_surface(int);
void pack_ny_surface(int);
void pack_nz_surface(int);
void pack_coordination(int);
void pack_viscosity(int);
void pack_pressure(int);
void pack_conductivity(int);
void pack_cv(int);
void pack_vx_shift(int);
void pack_vy_shift(int);
void pack_vz_shift(int);
};
} // namespace LAMMPS_NS
#endif
#endif

View File

@ -75,8 +75,8 @@ namespace RHEO_NS {
enum Status{
// Phase status
STATUS_FLUID = 1 << 0,
STATUS_REACTIVE = 1 << 1,
STATUS_SOLID = 1 << 2,
STATUS_SOLID = 1 << 1,
STATUS_REACTIVE = 1 << 2,
STATUS_FREEZING = 1 << 3,
// Surface status
STATUS_BULK = 1 << 4,

View File

@ -49,7 +49,7 @@ enum{TYPE,TYPE_FRACTION,TYPE_RATIO,TYPE_SUBSET,
DIPOLE,DIPOLE_RANDOM,SPIN_ATOM,SPIN_RANDOM,SPIN_ELECTRON,RADIUS_ELECTRON,
QUAT,QUAT_RANDOM,THETA,THETA_RANDOM,ANGMOM,OMEGA,TEMPERATURE,
DIAMETER,RADIUS_ATOM,DENSITY,VOLUME,IMAGE,BOND,ANGLE,DIHEDRAL,IMPROPER,
SPH_E,SPH_CV,SPH_RHO,EDPD_TEMP,EDPD_CV,CC,SMD_MASS_DENSITY,
RHEO_STATUS,SPH_E,SPH_CV,SPH_RHO,EDPD_TEMP,EDPD_CV,CC,SMD_MASS_DENSITY,
SMD_CONTACT_RADIUS,DPDTHETA,EPSILON,IVEC,DVEC,IARRAY,DARRAY};
#define BIG INT_MAX
@ -515,6 +515,24 @@ void Set::command(int narg, char **arg)
topology(IMPROPER);
iarg += 2;
} else if (strcmp(arg[iarg],"rheo/rho") == 0) {
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set rheo/rho", error);
if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1);
else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp);
if (!atom->rho_flag)
error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style());
set(SPH_RHO);
iarg += 2;
} else if (strcmp(arg[iarg],"rheo/status") == 0) {
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set rheo/status", error);
if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1);
else ivalue = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
if (!atom->status_flag)
error->all(FLERR,"Cannot set attribute {} for atom style {}", arg[iarg], atom->get_style());
set(RHEO_STATUS);
iarg += 2;
} else if (strcmp(arg[iarg],"sph/e") == 0) {
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "set sph/e", error);
if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1);
@ -881,6 +899,13 @@ void Set::set(int keyword)
if (dvalue <= 0.0) error->one(FLERR,"Invalid volume in set command");
atom->vfrac[i] = dvalue;
}
else if (keyword == RHEO_STATUS) {
if (ivalue != 0 && ivalue !=2)
error->one(FLERR,"Invalid value {} in set command for rheo/status", ivalue);
atom->status[i] = ivalue;
}
else if (keyword == SPH_E) atom->esph[i] = dvalue;
else if (keyword == SPH_CV) atom->cv[i] = dvalue;
else if (keyword == SPH_RHO) atom->rho[i] = dvalue;