ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Adding base fix and atom style

Browse Source
This commit is contained in:
jtclemm
2023-02-15 16:07:10 -07:00
parent 4cb29ce413
commit 348296e638
6 changed files with 661 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

135
src/RHEO/atom_vec_rheo.cpp Normal file
View File

@ -0,0 +1,135 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#include "atom_vec_rheo.h"
#include "atom.h"
#include <cstring>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
AtomVecRHEO::AtomVecRHEO(LAMMPS *lmp) : AtomVec(lmp)
{
molecular = Atom::ATOMIC;
mass_type = PER_TYPE;
forceclearflag = 1;
atom->status_flag = 1;
atom->rho_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 = {"status", "rho", "drho"};
fields_copy = {"status", "rho", "drho"};
fields_comm = {"status", "rho"};
fields_comm_vel = {"status", "rho"};
fields_reverse = {"drho"};
fields_border = {"status", "rho"};
fields_border_vel = {"status", "rho"};
fields_exchange = {"status", "rho"};
fields_restart = {"status", "rho"};
fields_create = {"status", "rho", "drho"};
fields_data_atom = {"id", "type", "status", "rho", "x"};
fields_data_vel = {"id", "v"};
setup_fields();
}
/* ----------------------------------------------------------------------
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 AtomVecRHEO::grow_pointers()
{
status = atom->status;
rho = atom->rho;
drho = atom->drho;
}
/* ----------------------------------------------------------------------
clear extra forces starting at atom N
nbytes = # of bytes to clear for a per-atom vector
------------------------------------------------------------------------- */
void AtomVecRHEO::force_clear(int n, size_t nbytes)
{
memset(&drho[n], 0, nbytes);
}
/* ----------------------------------------------------------------------
modify what AtomVec::data_atom() just unpacked
or initialize other atom quantities
------------------------------------------------------------------------- */
void AtomVecRHEO::data_atom_post(int ilocal)
{
drho[ilocal] = 0.0;
}
/* ----------------------------------------------------------------------
assign an index to named atom property and return index
return -1 if name is unknown to this atom style
------------------------------------------------------------------------- */
int AtomVecRHEO::property_atom(const std::string &name)
{
if (name == "status") return 0;
if (name == "rho") return 1;
if (name == "drho") return 2;
return -1;
}
/* ----------------------------------------------------------------------
pack per-atom data into buf for ComputePropertyAtom
index maps to data specific to this atom style
------------------------------------------------------------------------- */
void AtomVecRHEO::pack_property_atom(int index, double *buf, int nvalues, int groupbit)
{
int *mask = atom->mask;
int nlocal = atom->nlocal;
int n = 0;
if (index == 0) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = status[i];
else
buf[n] = 0.0;
n += nvalues;
}
} if else (index == 1) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = rho[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 2) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = drho[i];
else
buf[n] = 0.0;
n += nvalues;
}
}
}

45
src/RHEO/atom_vec_rheo.h Normal file
View File

@ -0,0 +1,45 @@
/* -*- 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 ATOM_CLASS
// clang-format off
AtomStyle(rheo,AtomVecRHEO);
// clang-format on
#else
#ifndef LMP_ATOM_VEC_RHEO_H
#define LMP_ATOM_VEC_RHEO_H
#include "atom_vec.h"
namespace LAMMPS_NS {
class AtomVecRHEO : virtual public AtomVec {
public:
AtomVecRHEO(class LAMMPS *);
void grow_pointers() override;
void force_clear(int, size_t) override;
void data_atom_post(int) override;
int property_atom(const std::string &) override;
void pack_property_atom(int, double *, int, int) override;
private:
int *status;
double *rho, *drho;
};
} // namespace LAMMPS_NS
#endif
#endif

365
src/RHEO/fix_rheo.cpp Normal file
View File

@ -0,0 +1,365 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, 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_rheo.h"
#include "atom.h"
#include "compute_rheo_grad.h"
#include "compute_rheo_interface.h"
#include "compute_rheo_kernel.h"
#include "compute_rheo_rhosum.h"
#include "compute_rheo_vshift.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "modify.h"
#include "update.h"
#include "utils.h"
using namespace LAMMPS_NS;
using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), compute_grad(nullptr), compute_kernel(nullptr),
compute_interface(nullptr), compute_rhosum(nullptr), compute_vshift(nullptr)
{
thermal_flag = 0;
rhosum_flag = 0;
shift_flag = 0;
solid_flag = 0;
rho0 = 1.0;
csq = 1.0;
if (atom->rho_flag != 1)
error->all(FLERR,"fix rheo command requires atom_style with density");
if (atom->status_flag != 1)
error->all(FLERR,"fix rheo command requires atom_style with status");
if (narg < 5)
error->all(FLERR,"Insufficient arguments for fix rheo command");
cut = utils::numeric(FLERR,arg[3],false,lmp);
if (strcmp(arg[4],"Quintic") == 0) {
kernel_style = QUINTIC;
} else if (strcmp(arg[4],"CRK0") == 0) {
kernel_style = CRK0;
} else if (strcmp(arg[4],"CRK1") == 0) {
kernel_style = CRK1;
} else if (strcmp(arg[4],"CRK2") == 0) {
kernel_style = CRK2;
} else error->all(FLERR,"Unknown kernel style {} in fix rheo", arg[4]);
zmin_kernel = utils::numeric(FLERR,arg[5],false,lmp);
int iarg = 6;
while (iarg < narg){
if (strcmp(arg[iarg],"shift") == 0) {
shift_flag = 1;
} else if (strcmp(arg[iarg],"thermal") == 0) {
thermal_flag = 1;
} else if (strcmp(arg[iarg],"rhosum") == 0) {
rhosum_flag = 1;
if(iarg + 1 >= narg) error->all(FLERR,"Illegal rhosum option in fix rheo");
zmin_rhosum = utils::inumeric(FLERR,arg[iarg + 1],false,lmp);
iarg += 1;
} else if (strcmp(arg[iarg],"rho0") == 0) {
if(iarg + 1 >= narg) error->all(FLERR,"Illegal rho0 option in fix rheo");
rho0 = utils::numeric(FLERR,arg[iarg + 1],false,lmp);
iarg += 1;
} else if (strcmp(arg[iarg],"csq") == 0) {
if(iarg+1 >= narg) error->all(FLERR,"Illegal csq option in fix rheo");
csq = utils::numeric(FLERR,arg[iarg + 1],false,lmp);
iarg += 1;
} else {
error->all(FLERR, "Illegal fix rheo command: {}", arg[iarg]);
}
iarg += 1;
}
time_integrate = 1;
thermal_fix_defined = 0;
viscosity_fix_defined = 0;
pressure_fix_defined = 0;
}
/* ---------------------------------------------------------------------- */
FixRHEO::~FixRHEO()
{
if (compute_kernel) modify->delete_compute("rheo_kernel");
if (compute_grad) modify->delete_compute("rheo_grad");
if (compute_interface) modify->delete_compute("rheo_interface");
if (compute_rhosum) modify->delete_compute("rheo_rhosum");
if (compute_vshift) modify->delete_compute("rheo_vshift");
}
/* ---------------------------------------------------------------------- */
void FixRHEO::post_constructor()
{
compute_kernel = dynamic_cast<ComputeRHEOKernel *>(modify->add_compute(fmt::format("rheo_kernel all rheo/kernel {} {} {}", kernel_style, zmin_kernel, cut)));
if (thermal_flag)
compute_grad = dynamic_cast<ComputeRHEOGrad *>(modify->add_compute(fmt::format("rheo_grad all rheo/grad {} velocity rho viscosity temprature", cut)));
else
compute_grad = dynamic_cast<ComputeRHEOGrad *>(modify->add_compute(fmt::format("rheo_grad all rheo/grad {} velocity rho viscosity", cut)));
compute_interface = dynamic_cast<ComputeRHEOInterface *>(modify->add_compute(fmt::format("rheo_interface all rheo/interface {}", cut)));
if (rhosum_flag)
compute_rhosum = dynamic_cast<ComputeRHEORhoSum *>(modify->add_compute(fmt::format("rheo_rhosum all rheo/rho/sum {} {}", cut, zmin_rhosum)));
if (shift_flag)
compute_vshift = dynamic_cast<ComputeRHEOVShift *>(modify->add_compute(fmt::format("rheo_vshift all rheo/vshift {}", cut)));
}
/* ---------------------------------------------------------------------- */
int FixRHEO::setmask()
{
int mask = 0;
mask |= INITIAL_INTEGRATE;
mask |= FINAL_INTEGRATE;
mask |= PRE_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixRHEO::init()
{
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
}
/* ---------------------------------------------------------------------- */
void FixRHEO::setup_pre_force(int /*vflag*/)
{
// Check to confirm all accessory fixes are defined
if(!thermal_fix_defined && thermal_flag)
error->all(FLERR, "Missing fix rheo/thermal");
if (!viscosity_fix_defined)
error->all(FLERR, "Missing fix rheo/viscosity");
if (!pressure_fix_defined)
error->all(FLERR, "Missing fix rheo/pressure");
// Reset to zero for next run
thermal_fix_defined = 0;
viscosity_fix_defined = 0;
pressure_fix_defined = 0;
pre_force(0);
}
/* ---------------------------------------------------------------------- */
void FixRHEO::initial_integrate(int /*vflag*/)
{
// update v and x and rho of atoms in group
int i, a, b;
double dtfm, divu;
int dim = domain->dimension;
int *status = atom->status;
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *rho = atom->rho;
double *drho = atom->drho;
double *mass = atom->mass;
double *rmass = atom->rmass;
int rmass_flag = atom->rmass_flag;
double **gradr = compute_grad->gradr;
double **gradv = compute_grad->gradv;
double **vshift = compute_vshift->array_atom;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup)
nlocal = atom->nfirst;
//Density Half-step
for (i = 0; i < nlocal; i++) {
if (status[i] & STATUS_NO_FORCE) continue;
if (mask[i] & groupbit) {
if (rmass_flag) {
dtfm = dtf / rmass[i];
} else {
dtfm = dtf / mass[type[i]];
}
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
}
}
// Update gradients and interpolate solid properties
compute_grad->forward_fields(); // also forwards v and rho for chi
compute_interface->store_forces(); // Need to save, wiped in exchange
compute_interface->compute_peratom();
compute_grad->compute_peratom();
// Position half-step
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
for (a = 0; a < dim; a++) {
x[i][a] += dtv * v[i][a];
}
}
}
// Update density using div(u)
if (!rhosum_flag) {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (status[i] & STATUS_NO_FORCE) continue;
if (!(status[i] & STATUS_FLUID)) continue;
divu = 0;
for (a = 0; a < dim; a++) {
divu += gradv[i][a * (1 + dim)];
}
rho[i] += dtf * (drho[i] - rho[i] * divu);
}
}
}
// Shifting atoms
if (shift_flag) {
compute_vshift->correct_surfaces();
for (i = 0; i < nlocal; i++) {
if (!(status[i] & STATUS_SHIFT)) continue;
if (mask[i] & groupbit) {
for (a = 0; a < dim; a++) {
x[i][a] += dtv * vshift[i][a];
for (b = 0; b < dim; b++) {
v[i][a] += dtv * vshift[i][b] * gradv[i][a * dim + b];
}
}
if (!rhosum_flag) {
for (a = 0; a < dim; a++) {
rho[i] += dtv * vshift[i][a] * gradr[i][a];
}
}
}
}
}
}
/* ---------------------------------------------------------------------- */
void FixRHEO::pre_force(int /*vflag*/)
{
if (rhosum_flag)
compute_rhosum->compute_peratom();
compute_grad->forward_fields(); // also forwards v and rho for chi
compute_kernel->compute_peratom();
compute_interface->compute_peratom();
compute_grad->compute_peratom();
compute_grad->forward_gradients();
if (shift_flag)
compute_vshift->compute_peratom();
// Remove extra shifting/no force options options
int *status = atom->status;
int nall = atom->nlocal + atom->nghost;
for (int i = 0; i < nall; i++) {
if (mask[i] & groupbit) {
status[i] &= ~STATUS_NO_FORCE;
if (status[i] & STATUS_FLUID)
status[i] &= ~STATUS_SHIFT;
}
}
}
/* ---------------------------------------------------------------------- */
void FixRHEO::final_integrate() {
int *status = atom->status;
double **gradv = compute_grad->gradv;
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *rho = atom->rho;
double *drho = atom->drho;
int *type = atom->type;
int *mask = atom->mask;
double *mass = atom->mass;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup)
nlocal = atom->nfirst;
double dtfm, divu;
double *rmass = atom->rmass;
int rmass_flag = atom->rmass_flag;
int i, a;
int dim = domain->dimension;
// Update velocity
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (status[i] & STATUS_NO_FORCE) continue;
if (rmass_flag) {
dtfm = dtf / rmass[i];
} else {
dtfm = dtf / mass[type[i]];
}
for (a = 0; a < dim; a++) {
v[i][a] += dtfm * f[i][a];
}
}
}
// Update density using divu
if (!rhosum_flag) {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (status[i] & STATUS_NO_FORCE) continue;
if (!(status[i] & STATUS_FLUID)) continue;
divu = 0;
for (a = 0; a < dim; a++) {
divu += gradv[i][a * (1 + dim)];
}
rho[i] += dtf * (drho[i] - rho[i] * divu);
}
}
}
}
/* ---------------------------------------------------------------------- */
void FixRHEO::reset_dt()
{
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
}

91
src/RHEO/fix_rheo.h Normal file
View File

@ -0,0 +1,91 @@
/* -*- 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 FIX_CLASS
// clang-format off
FixStyle(rheo,FixRHEO)
// clang-format on
#else
#ifndef LMP_FIX_RHEO_H
#define LMP_FIX_RHEO_H
#include "fix.h"
namespace LAMMPS_NS {
class FixRHEO : public Fix {
public:
FixRHEO(class LAMMPS *, int, char **);
virtual ~FixRHEO();
int setmask();
virtual void post_constructor();
virtual void init();
virtual void setup_pre_force(int);
virtual void pre_force(int);
virtual void initial_integrate(int);
virtual void final_integrate();
void reset_dt();
int kernel_style;
int thermal_flag;
int rhosum_flag;
int shift_flag;
int solid_flag;
int thermal_fix_defined;
int viscosity_fix_defined;
int pressure_fix_defined;
int *status, *surface;
double *conductivity, *viscosity, *pressure;
double **f_pressure;
class ComputeRHEOGrad *compute_grad;
class ComputeRHEOKernel *compute_kernel;
class ComputeRHEOInterface *compute_interface;
class ComputeRHEORhoSum *compute_rhosum;
class ComputeRHEOVShift *compute_vshift;
enum {QUINTIC, CRK0, CRK1, CRK2};
enum {LINEAR, CUBIC, TAITWATER};
enum {
// Phase status
STATUS_FLUID = 1 << 0,
STATUS_REACTIVE = 1 << 1,
STATUS_SOLID = 1 << 2,
STATUS_FREEZING = 1 << 3
// Temporary status options - reset in preforce
STATUS_SHIFT = 1 << 4,
STATUS_NO_FORCE = 1 << 5,
// Surface status
STATUS_BULK = 1 << 6,
STATUS_LAYER = 1 << 7,
STATUS_SURFACE = 1 << 8,
STATUS_SPLASH = 1 << 9,
};
protected:
double cut, rho0, csq;
int zmin_kernel, rhosum_zmin;
double dtv, dtf;
};
} // namespace LAMMPS_NS
#endif
#endif

20
src/atom.cpp
View File

@ -195,6 +195,10 @@ Atom::Atom(LAMMPS *_lmp) : Pointers(_lmp)
eff_plastic_strain_rate = nullptr;
damage = nullptr;
// RHEO package
status = nullptr;
// SPH package
rho = drho = esph = desph = cv = nullptr;
@ -521,6 +525,10 @@ void Atom::peratom_create()
add_peratom("cc",&cc,DOUBLE,1);
add_peratom("cc_flux",&cc_flux,DOUBLE,1,1); // set per-thread flag
// RHEO package
add_peratom("status",&status,INT,0);
// SPH package
add_peratom("rho",&rho,DOUBLE,0);
@ -625,6 +633,7 @@ void Atom::set_atomflag_defaults()
rmass_flag = radius_flag = omega_flag = torque_flag = angmom_flag = 0;
vfrac_flag = spin_flag = eradius_flag = ervel_flag = erforce_flag = 0;
cs_flag = csforce_flag = vforce_flag = ervelforce_flag = etag_flag = 0;
status_flag = 0;
rho_flag = esph_flag = cv_flag = vest_flag = 0;
dpd_flag = edpd_flag = tdpd_flag = 0;
sp_flag = 0;
@ -2911,7 +2920,12 @@ void *Atom::extract(const char *name)
if (strcmp(name,"vforce") == 0) return (void *) vforce;
if (strcmp(name,"etag") == 0) return (void *) etag;
// RHEO package
if (strcmp(name,"status") == 0) return (void *) status;
// SPH package
if (strcmp(name,"rho") == 0) return (void *) rho;
if (strcmp(name,"drho") == 0) return (void *) drho;
if (strcmp(name,"esph") == 0) return (void *) esph;
@ -3030,6 +3044,12 @@ int Atom::extract_datatype(const char *name)
if (strcmp(name,"vforce") == 0) return LAMMPS_DOUBLE_2D;
if (strcmp(name,"etag") == 0) return LAMMPS_INT;
// RHEO package
if (strcmp(name,"status") == 0) return LAMMPS_INT;
// SPH package
if (strcmp(name,"rho") == 0) return LAMMPS_DOUBLE;
if (strcmp(name,"drho") == 0) return LAMMPS_DOUBLE;
if (strcmp(name,"esph") == 0) return LAMMPS_DOUBLE;

5
src/atom.h
View File

@ -154,6 +154,10 @@ class Atom : protected Pointers {
double *eff_plastic_strain_rate;
double *damage;
// RHEO package
int *status;
// SPH package
double *rho, *drho, *esph, *desph, *cv;
@ -188,6 +192,7 @@ class Atom : protected Pointers {
int rmass_flag, radius_flag, omega_flag, torque_flag, angmom_flag, quat_flag;
int vfrac_flag, spin_flag, eradius_flag, ervel_flag, erforce_flag;
int cs_flag, csforce_flag, vforce_flag, ervelforce_flag, etag_flag;
int status_flag;
int rho_flag, esph_flag, cv_flag, vest_flag;
int dpd_flag, edpd_flag, tdpd_flag;
int mesont_flag;