From c6d59fc5267bede32772403f683c136e36d8528a Mon Sep 17 00:00:00 2001 From: jtclemm Date: Sun, 12 Jun 2022 10:09:26 -0600 Subject: [PATCH] adding documentation and integration fix --- src/.gitignore | 6 +- src/GRANULAR/atom_vec_sphere_temperature.cpp | 160 +++++++++++++++++++ src/GRANULAR/atom_vec_sphere_temperature.h | 51 ++++++ src/GRANULAR/fix_temp_integrate.cpp | 142 ++++++++++++++++ src/GRANULAR/fix_temp_integrate.h | 50 ++++++ src/atom.cpp | 18 +++ src/atom.h | 2 + src/set.cpp | 19 ++- 8 files changed, 445 insertions(+), 3 deletions(-) create mode 100644 src/GRANULAR/atom_vec_sphere_temperature.cpp create mode 100644 src/GRANULAR/atom_vec_sphere_temperature.h create mode 100644 src/GRANULAR/fix_temp_integrate.cpp create mode 100644 src/GRANULAR/fix_temp_integrate.h diff --git a/src/.gitignore b/src/.gitignore index 6657256e8f..7a6bc1beb4 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -421,14 +421,14 @@ /atom_vec_full.h /atom_vec_full_hars.cpp /atom_vec_full_hars.h -/atom_vec_granular.cpp -/atom_vec_granular.h /atom_vec_molecular.cpp /atom_vec_molecular.h /atom_vec_oxdna.cpp /atom_vec_oxdna.h /atom_vec_peri.cpp /atom_vec_peri.h +/atom_vec_sphere_temperature.cpp +/atom_vec_sphere_temperature.h /atom_vec_template.cpp /atom_vec_template.h /body_nparticle.cpp @@ -1510,6 +1510,8 @@ /fix_smd_wall_surface.h /fix_srp.cpp /fix_srp.h +/fix_temp_integrate.cpp +/fix_temp_integrate.h /fix_tfmc.cpp /fix_tfmc.h /fix_ttm.cpp diff --git a/src/GRANULAR/atom_vec_sphere_temperature.cpp b/src/GRANULAR/atom_vec_sphere_temperature.cpp new file mode 100644 index 0000000000..5b62a2c279 --- /dev/null +++ b/src/GRANULAR/atom_vec_sphere_temperature.cpp @@ -0,0 +1,160 @@ +/* ---------------------------------------------------------------------- + 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 "atom_vec_sphere_temperature.h" + +#include "atom.h" +#include "error.h" +#include "fix.h" +#include "fix_adapt.h" +#include "math_const.h" +#include "modify.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +/* ---------------------------------------------------------------------- */ + +AtomVecSphereTemperature::AtomVecSphereTemperature(LAMMPS *lmp) : AtomVec(lmp) +{ + mass_type = PER_ATOM; + molecular = Atom::ATOMIC; + + atom->sphere_flag = 1; + atom->radius_flag = atom->rmass_flag = atom->omega_flag = atom->torque_flag = 1; + atom->temperature_flag = atom->heatflux_flag = 1; + + // strings with peratom variables to include in each AtomVec method + // strings cannot contain fields in corresponding AtomVec default strings + // order of fields in a string does not matter + // except: fields_data_atom & fields_data_vel must match data file + + fields_grow = {"radius", "rmass", "omega", "torque", "temperature", "heatflux"}; + fields_copy = {"radius", "rmass", "omega", "temperature"}; + fields_comm_vel = {"omega", "temperature"}; + fields_reverse = {"torque", "heatflux"}; + fields_border = {"radius", "rmass", "temperature"}; + fields_border_vel = {"radius", "rmass", "omega", "temperature"}; + fields_exchange = {"radius", "rmass", "omega", "temperature"}; + fields_restart = {"radius", "rmass", "omega", "temperature"}; + fields_create = {"radius", "rmass", "omega", "temperature"}; + fields_data_atom = {"id", "type", "radius", "rmass", "x", "temperature"}; + fields_data_vel = {"id", "v", "omega"}; +} + +/* ---------------------------------------------------------------------- + process sub-style args + optional arg = 0/1 for static/dynamic particle radii +------------------------------------------------------------------------- */ + +void AtomVecSphereTemperature::process_args(int narg, char **arg) +{ + if (narg != 0 && narg != 1) error->all(FLERR, "Illegal atom_style sphere command"); + + radvary = 0; + if (narg == 1) { + radvary = utils::numeric(FLERR, arg[0], true, lmp); + if (radvary < 0 || radvary > 1) error->all(FLERR, "Illegal atom_style sphere command"); + } + + // dynamic particle radius and mass must be communicated every step + + if (radvary) { + fields_comm = {"radius", "rmass"}; + fields_comm_vel = {"radius", "rmass", "omega"}; + } + + // delay setting up of fields until now + + setup_fields(); +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecSphereTemperature::init() +{ + AtomVec::init(); + + // check if optional radvary setting should have been set to 1 + + for (int i = 0; i < modify->nfix; i++) + if (strcmp(modify->fix[i]->style, "adapt") == 0) { + auto fix = dynamic_cast(modify->fix[i]); + if (fix->diamflag && radvary == 0) + error->all(FLERR, "Fix adapt changes particle radii but atom_style sphere is not dynamic"); + } +} + +/* ---------------------------------------------------------------------- + set local copies of all grow ptrs used by this class, except defaults + needed in replicate when 2 atom classes exist and it calls pack_restart() +------------------------------------------------------------------------- */ + +void AtomVecSphereTemperature::grow_pointers() +{ + radius = atom->radius; + rmass = atom->rmass; + omega = atom->omega; +} + +/* ---------------------------------------------------------------------- + initialize non-zero atom quantities +------------------------------------------------------------------------- */ + +void AtomVecSphereTemperature::create_atom_post(int ilocal) +{ + radius[ilocal] = 0.5; + rmass[ilocal] = 4.0 * MY_PI / 3.0 * 0.5 * 0.5 * 0.5; +} + +/* ---------------------------------------------------------------------- + modify what AtomVec::data_atom() just unpacked + or initialize other atom quantities +------------------------------------------------------------------------- */ + +void AtomVecSphereTemperature::data_atom_post(int ilocal) +{ + radius_one = 0.5 * atom->radius[ilocal]; + radius[ilocal] = radius_one; + if (radius_one > 0.0) rmass[ilocal] *= 4.0 * MY_PI / 3.0 * radius_one * radius_one * radius_one; + + if (rmass[ilocal] <= 0.0) error->one(FLERR, "Invalid density in Atoms section of data file"); + + omega[ilocal][0] = 0.0; + omega[ilocal][1] = 0.0; + omega[ilocal][2] = 0.0; +} + +/* ---------------------------------------------------------------------- + modify values for AtomVec::pack_data() to pack +------------------------------------------------------------------------- */ + +void AtomVecSphereTemperature::pack_data_pre(int ilocal) +{ + radius_one = radius[ilocal]; + rmass_one = rmass[ilocal]; + + radius[ilocal] *= 2.0; + if (radius_one != 0.0) + rmass[ilocal] = rmass_one / (4.0 * MY_PI / 3.0 * radius_one * radius_one * radius_one); +} + +/* ---------------------------------------------------------------------- + unmodify values packed by AtomVec::pack_data() +------------------------------------------------------------------------- */ + +void AtomVecSphereTemperature::pack_data_post(int ilocal) +{ + radius[ilocal] = radius_one; + rmass[ilocal] = rmass_one; +} diff --git a/src/GRANULAR/atom_vec_sphere_temperature.h b/src/GRANULAR/atom_vec_sphere_temperature.h new file mode 100644 index 0000000000..6b1ca6f390 --- /dev/null +++ b/src/GRANULAR/atom_vec_sphere_temperature.h @@ -0,0 +1,51 @@ +/* -*- 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 ATOM_CLASS +// clang-format off +AtomStyle(sphere,AtomVecSphereTemperature); +// clang-format on +#else + +#ifndef LMP_ATOM_VEC_SPHERE_TEMPERATURE_H +#define LMP_ATOM_VEC_SPHERE_TEMPERATURE_H + +#include "atom_vec.h" + +namespace LAMMPS_NS { + +class AtomVecSphereTemperature : public AtomVec { + public: + AtomVecSphereTemperature(class LAMMPS *); + void process_args(int, char **) override; + void init() override; + + void grow_pointers() override; + void create_atom_post(int) override; + void data_atom_post(int) override; + void pack_data_pre(int) override; + void pack_data_post(int) override; + + private: + double *radius, *rmass; + double **omega; + double *temperature, *heatflux; + + int radvary; + double radius_one, rmass_one; +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/GRANULAR/fix_temp_integrate.cpp b/src/GRANULAR/fix_temp_integrate.cpp new file mode 100644 index 0000000000..bb351d7db4 --- /dev/null +++ b/src/GRANULAR/fix_temp_integrate.cpp @@ -0,0 +1,142 @@ +// 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_temp_integrate.h" + +#include "atom.h" +#include "error.h" +#include "force.h" +#include "memory.h" +#include "respa.h" +#include "update.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +enum {NONE, CONSTANT, TYPE}; + +/* ---------------------------------------------------------------------- */ + +FixTempIntegrate::FixTempIntegrate(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + if (narg < 4) error->all(FLERR,"Illegal fix command"); + + cp_style = NONE; + + int ntypes = atom->ntypes; + int iarg = 3; + while (iarg < narg) { + if (strcmp(arg[iarg+1],"constant") == 0) { + if (iarg+2 >= narg) error->all(FLERR,"Illegal fix command"); + cp_style = CONSTANT; + cp = utils::numeric(FLERR,arg[iarg+2],false,lmp); + if (cp < 0.0) error->all(FLERR,"Illegal fix command"); + iarg += 2; + } else if (strcmp(arg[iarg+1],"type") == 0) { + if (iarg+1+ntypes >= narg) error->all(FLERR,"Illegal fix command"); + cp_style = TYPE; + memory->create(cp_type,ntypes+1,"fix/temp/integrate:cp_type"); + for (int i = 1; i <= ntypes; i++) { + cp_type[i] = utils::numeric(FLERR,arg[iarg+1+i],false,lmp); + if (cp_type[i] < 0.0) error->all(FLERR,"Illegal fix command"); + } + iarg += 1+ntypes; + } else { + error->all(FLERR,"Illegal fix command"); + } + iarg += 1; + } + + if (cp_style == NONE) + error->all(FLERR, "Must specify specific heat in fix temp/integrate"); + dynamic_group_allow = 1; + time_integrate = 1; +} + +/* ---------------------------------------------------------------------- */ + +int FixTempIntegrate::setmask() +{ + int mask = 0; + mask |= FINAL_INTEGRATE; + mask |= FINAL_INTEGRATE_RESPA; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixTempIntegrate::init() +{ + dt = update->dt; + + if (!atom->temperature_flag) + error->all(FLERR,"Fix temp/integrate requires atom style with temperature property"); + if (!atom->heatflux_flag) + error->all(FLERR,"Fix temp/integrate requires atom style with heatflux property"); +} + +/* ---------------------------------------------------------------------- */ + +void FixTempIntegrate::final_integrate() +{ + // update temperature of atoms in group + + double *temperature = atom->temperature; + double *heatflux = atom->heatflux; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + if (rmass) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + temperature[i] += dt * heatflux[i] / (calc_cp(i) * rmass[i]); + } + } else { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + temperature[i] += dt * heatflux[i] / (calc_cp(i) * mass[type[i]]); + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixTempIntegrate::final_integrate_respa(int ilevel, int /*iloop*/) +{ + dt = update->dt; + final_integrate(); +} + +/* ---------------------------------------------------------------------- */ + +void FixTempIntegrate::reset_dt() +{ + dt = update->dt; +} + +/* ---------------------------------------------------------------------- */ + +double FixTempIntegrate::calc_cp(int i) +{ + if (cp_style == TYPE) { + return cp_type[atom->type[i]]; + } else { + return cp; + } +} diff --git a/src/GRANULAR/fix_temp_integrate.h b/src/GRANULAR/fix_temp_integrate.h new file mode 100644 index 0000000000..738cb28a6f --- /dev/null +++ b/src/GRANULAR/fix_temp_integrate.h @@ -0,0 +1,50 @@ +/* -*- 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(temp/integrate,FixTempIntegrate); +// clang-format on +#else + +#ifndef LMP_FIX_TEMP_INTEGRATE_H +#define LMP_FIX_TEMP_INTEGRATE_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixTempIntegrate : public Fix { + public: + FixTempIntegrate(class LAMMPS *, int, char **); + + int setmask() override; + void init() override; + void final_integrate() override; + void final_integrate_respa(int, int) override; + void reset_dt() override; + + protected: + double dt; + double cp, *cp_type; + int cp_style; + + int mass_require; + + double calc_cp(int); +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/atom.cpp b/src/atom.cpp index 8b3ab8c75d..ce66742779 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -123,6 +123,8 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) radius = rmass = nullptr; ellipsoid = line = tri = body = nullptr; quat = nullptr; + temperature = nullptr; + heatflux = nullptr; // molecular systems @@ -410,6 +412,9 @@ void Atom::peratom_create() add_peratom("tri",&tri,INT,0); add_peratom("body",&body,INT,0); + add_peratom("temperature",&temperature,DOUBLE,0); + add_peratom("heatflux",&heatflux,DOUBLE,0); + // BPM package add_peratom("quat",&quat,DOUBLE,4); @@ -612,6 +617,7 @@ void Atom::set_atomflag_defaults() molecule_flag = molindex_flag = molatom_flag = 0; q_flag = mu_flag = 0; rmass_flag = radius_flag = omega_flag = torque_flag = angmom_flag = 0; + temperature_flag = heatflux_flag = 0; vfrac_flag = spin_flag = eradius_flag = ervel_flag = erforce_flag = 0; cs_flag = csforce_flag = vforce_flag = ervelforce_flag = etag_flag = 0; rho_flag = esph_flag = cv_flag = vest_flag = 0; @@ -2602,6 +2608,14 @@ length of the data area, and a short description. - double - 4 - four quaternion components of the particles + * - temperature + - double + - 1 + - temperature of the particles + * - heatflux + - double + - 1 + - heatflux of the particles * - i_name - int - 1 @@ -2658,6 +2672,8 @@ void *Atom::extract(const char *name) if (strcmp(name,"tri") == 0) return (void *) tri; if (strcmp(name,"body") == 0) return (void *) body; if (strcmp(name,"quat") == 0) return (void *) quat; + if (strcmp(name,"temperature") == 0) return (void *) temperature; + if (strcmp(name,"heatflux") == 0) return (void *) heatflux; if (strcmp(name,"vfrac") == 0) return (void *) vfrac; if (strcmp(name,"s0") == 0) return (void *) s0; @@ -2781,6 +2797,8 @@ int Atom::extract_datatype(const char *name) if (strcmp(name,"tri") == 0) return LAMMPS_INT; if (strcmp(name,"body") == 0) return LAMMPS_INT; if (strcmp(name,"quat") == 0) return LAMMPS_DOUBLE_2D; + if (strcmp(name,"temperature") == 0) return LAMMPS_DOUBLE; + if (strcmp(name,"heatflux") == 0) return LAMMPS_DOUBLE; if (strcmp(name,"vfrac") == 0) return LAMMPS_DOUBLE; if (strcmp(name,"s0") == 0) return LAMMPS_DOUBLE; diff --git a/src/atom.h b/src/atom.h index 499e5f0d0f..7ac459317c 100644 --- a/src/atom.h +++ b/src/atom.h @@ -80,6 +80,7 @@ class Atom : protected Pointers { double **omega, **angmom, **torque; int *ellipsoid, *line, *tri, *body; double **quat; + double *temperature, *heatflux; // molecular systems @@ -182,6 +183,7 @@ class Atom : protected Pointers { int molecule_flag, molindex_flag, molatom_flag; int q_flag, mu_flag; int rmass_flag, radius_flag, omega_flag, torque_flag, angmom_flag, quat_flag; + int temperature_flag, heatflux_flag; int vfrac_flag, spin_flag, eradius_flag, ervel_flag, erforce_flag; int cs_flag, csforce_flag, vforce_flag, ervelforce_flag, etag_flag; int rho_flag, esph_flag, cv_flag, vest_flag; diff --git a/src/set.cpp b/src/set.cpp index 2843281d78..7c759873a2 100644 --- a/src/set.cpp +++ b/src/set.cpp @@ -47,7 +47,7 @@ enum{ATOM_SELECT,MOL_SELECT,TYPE_SELECT,GROUP_SELECT,REGION_SELECT}; enum{TYPE,TYPE_FRACTION,TYPE_RATIO,TYPE_SUBSET, MOLECULE,X,Y,Z,VX,VY,VZ,CHARGE,MASS,SHAPE,LENGTH,TRI, DIPOLE,DIPOLE_RANDOM,SPIN,SPIN_RANDOM,QUAT,QUAT_RANDOM, - THETA,THETA_RANDOM,ANGMOM,OMEGA, + THETA,THETA_RANDOM,ANGMOM,OMEGA,TEMPERATURE, DIAMETER,DENSITY,VOLUME,IMAGE,BOND,ANGLE,DIHEDRAL,IMPROPER, SPH_E,SPH_CV,SPH_RHO,EDPD_TEMP,EDPD_CV,CC,SMD_MASS_DENSITY, SMD_CONTACT_RADIUS,DPDTHETA,EPSILON,IVEC,DVEC,IARRAY,DARRAY}; @@ -389,6 +389,15 @@ void Set::command(int narg, char **arg) set(DENSITY); iarg += 2; + } else if (strcmp(arg[iarg],"temperature") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal set command"); + if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); + else dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (!atom->temperature_flag) + error->all(FLERR,"Cannot set this attribute for this atom style"); + set(TEMPERATURE); + iarg += 2; + } else if (strcmp(arg[iarg],"volume") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal set command"); if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1); @@ -768,6 +777,7 @@ void Set::set(int keyword) case SHAPE: case DIAMETER: case DENSITY: + case TEMPERATURE: case QUAT: case IMAGE: if (modify->check_rigid_list_overlap(select)) @@ -1008,6 +1018,13 @@ void Set::set(int keyword) atom->omega[i][2] = zvalue; } + // set temperature of particle + + else if (keyword == ANGMOM) { + if (dvalue < 0.0) error->one(FLERR,"Invalid temperature in set command"); + atom->temperature[i] = dvalue; + } + // reset any or all of 3 image flags else if (keyword == IMAGE) {