From 560bd90e1149bd880600b1edea291b7b2f1963cb Mon Sep 17 00:00:00 2001 From: jtclemm Date: Sun, 19 Feb 2023 22:08:02 -0700 Subject: [PATCH] Drafting viscosity fix --- src/RHEO/compute_rheo_grad.cpp | 11 +- src/RHEO/compute_rheo_grad.h | 16 +-- src/RHEO/fix_rheo.cpp | 9 ++ src/RHEO/fix_rheo.h | 3 +- src/RHEO/fix_rheo_viscosity.cpp | 198 ++++++++++++++++++++++++++++++++ src/RHEO/fix_rheo_viscosity.h | 49 ++++++++ 6 files changed, 272 insertions(+), 14 deletions(-) create mode 100644 src/RHEO/fix_rheo_viscosity.cpp create mode 100644 src/RHEO/fix_rheo_viscosity.h diff --git a/src/RHEO/compute_rheo_grad.cpp b/src/RHEO/compute_rheo_grad.cpp index 33dd3dc3bc..e87d39aa53 100644 --- a/src/RHEO/compute_rheo_grad.cpp +++ b/src/RHEO/compute_rheo_grad.cpp @@ -23,11 +23,8 @@ #include "force.h" #include "neighbor.h" #include "neigh_list.h" -#include "pair.h" -#include "memory.h" #include "modify.h" #include "update.h" -#include "utils.h" #include #include @@ -52,14 +49,13 @@ ComputeRHEOGrad::ComputeRHEOGrad(LAMMPS *lmp, int narg, char **arg) : else error->all(FLERR, "Illegal compute rheo/grad command, {}", arg[iarg]); } - dim = domain->dimension; - ncomm_grad = 0; ncomm_field = 0; comm_reverse = 0; std::string fix_cmd = "rheo_grad_property_atom all property/atom" + int dim = domain->dimension; if (velocity_flag) { ncomm_grad += dim * dim; ncomm_field += dim; @@ -165,6 +161,7 @@ void ComputeRHEOGrad::compute_peratom() int *type = atom->type; double *mass = atom->mass; int newton = force->newton; + int dim = domain->dimension; inum = list->inum; ilist = list->ilist; @@ -310,6 +307,7 @@ int ComputeRHEOGrad::pack_forward_comm(int n, int *list, double *buf, double *temperature = atom->temperature; double *eta = atom->viscosity; double **v = atom->v; + int dim = domain->dimension; m = 0; @@ -363,6 +361,7 @@ void ComputeRHEOGrad::unpack_forward_comm(int n, int first, double *buf) double * rho = atom->rho; double * temperature = atom->temperature; double ** v = atom->v; + int dim = domain->dimension; m = 0; last = first + n; @@ -404,6 +403,7 @@ void ComputeRHEOGrad::unpack_forward_comm(int n, int first, double *buf) int ComputeRHEOGrad::pack_reverse_comm(int n, int first, double *buf) { int i,k,m,last; + int dim = domain->dimension; m = 0; last = first + n; @@ -433,6 +433,7 @@ int ComputeRHEOGrad::pack_reverse_comm(int n, int first, double *buf) void ComputeRHEOGrad::unpack_reverse_comm(int n, int *list, double *buf) { int i,k,j,m; + int dim = domain->dimension; m = 0; for (i = 0; i < n; i++) { diff --git a/src/RHEO/compute_rheo_grad.h b/src/RHEO/compute_rheo_grad.h index 8b8e28fd98..35b50ca834 100644 --- a/src/RHEO/compute_rheo_grad.h +++ b/src/RHEO/compute_rheo_grad.h @@ -28,13 +28,13 @@ class ComputeRHEOGrad : public Compute { public: ComputeRHEOGrad(class LAMMPS *, int, char **); ~ComputeRHEOGrad(); - void init(); - void init_list(int, class NeighList *); - void compute_peratom(); - int pack_forward_comm(int, int *, double *, int, int *); - void unpack_forward_comm(int, int, double *); - int pack_reverse_comm(int, int, double *); - void unpack_reverse_comm(int, int *, double *); + void init() override; + void init_list(int, class NeighList *) override; + void compute_peratom() override; + int pack_forward_comm(int, int *, double *, int, int *) override; + void unpack_forward_comm(int, int, double *) override; + int pack_reverse_comm(int, int, double *) override; + void unpack_reverse_comm(int, int *, double *) override; void forward_gradients(); void forward_fields(); double **gradv; @@ -44,7 +44,7 @@ class ComputeRHEOGrad : public Compute { int stage; private: - int dim, comm_stage; + int comm_stage; int ncomm_grad, ncomm_field; double cut, cutsq, rho0; class NeighList *list; diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index 17108b4c91..b77a35dd95 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -123,6 +123,12 @@ void FixRHEO::post_constructor() if (shift_flag) compute_vshift = dynamic_cast(modify->add_compute(fmt::format("rheo_vshift all rheo/vshift {}", cut))); + + + //todo here + //allocate memory for viscosity, pressure, etc + //don't want to save to restart/datafiles (could disable fix store/state) + //but do want it available for dupm files } /* ---------------------------------------------------------------------- */ @@ -142,6 +148,9 @@ void FixRHEO::init() { dtv = update->dt; dtf = 0.5 * update->dt * force->ftm2v; + + if (modify->get_fix_by_style("rheo").size() > 1) + error->all(FLERR,"Can only specify one instance of fix rheo"); } /* ---------------------------------------------------------------------- */ diff --git a/src/RHEO/fix_rheo.h b/src/RHEO/fix_rheo.h index 8597eb2809..24be629efe 100644 --- a/src/RHEO/fix_rheo.h +++ b/src/RHEO/fix_rheo.h @@ -47,7 +47,8 @@ class FixRHEO : public Fix { int viscosity_fix_defined; int pressure_fix_defined; - int *status, *surface; + // Non-persistent per-atom arrays are initialized here + int *surface; double *conductivity, *viscosity, *pressure; double **f_pressure; diff --git a/src/RHEO/fix_rheo_viscosity.cpp b/src/RHEO/fix_rheo_viscosity.cpp new file mode 100644 index 0000000000..9f7d3cf23c --- /dev/null +++ b/src/RHEO/fix_rheo_viscosity.cpp @@ -0,0 +1,198 @@ +/* ---------------------------------------------------------------------- + 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_viscosity.h" + +#include "atom.h" +#include "comm.h" +#include "compute_rheo_grad.h" +#include "domain.h" +#include "error.h" +#include "fix_rheo.h" +#include "memory.h" +#include "modify.h" +#include "update.h" + +#include + +using namespace LAMMPS_NS; +using namespace FixConst; +enum {NONE, CONSTANT, TYPE, POWER}; + +/* ---------------------------------------------------------------------- */ + +FixRHEOViscosity::FixRHEOViscosity(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg), eta_type(nullptr) +{ + if (narg < 4) error->all(FLERR,"Illegal fix command"); + + viscosity_style = NONE; + + comm_forward = 1; + + int ntypes = atom->ntypes; + int iarg = 3; + if (strcmp(arg[iarg],"constant") == 0) { + if (iarg+1 >= narg) error->all(FLERR,"Insufficient arguments for fix rheo/viscosity"); + viscosity_style = CONSTANT; + eta = utils::numeric(FLERR,arg[iarg + 1],false,lmp); + if (eta < 0.0) error->all(FLERR,"The viscosity must be positive in fix rheo/viscosity"); + iarg += 1; + } else if (strcmp(arg[iarg],"type") == 0) { + if(iarg+ntypes >= narg) error->all(FLERR,"Insufficient arguments for fix rheo/viscosity"); + viscosity_style = TYPE; + memory->create(eta_type, ntypes + 1, "rheo_thermal:eta_type"); + for (int i = 1; i <= ntypes; i++) { + eta_type[i] = utils::numeric(FLERR,arg[iarg + 1 + i], false, lmp); + if (eta_type[i] < 0.0) error->all(FLERR,"The viscosity must be positive in fix rheo/viscosity"); + } + iarg += ntypes; + } else if (strcmp(arg[iarg],"power") == 0) { + if (iarg+4 >= narg) error->all(FLERR,"Insufficient arguments for fix rheo/viscosity"); + viscosity_style = POWER; + eta = utils::numeric(FLERR,arg[iarg + 1],false,lmp); + gd0 = utils::numeric(FLERR,arg[iarg + 2],false,lmp); + K = utils::numeric(FLERR,arg[iarg + 3],false,lmp); + npow = utils::numeric(FLERR,arg[iarg + 4],false,lmp); + tau0 = eta * gd0 - K * pow(gd0, npow); + if (eta < 0.0) error->all(FLERR,"The viscosity must be positive in fix rheo/viscosity"); + iarg += 5; + } else { + error->all(FLERR,"Illegal fix command, {}", arg[iarg]); + } + + if (viscosity_style == NONE) + error->all(FLERR,"Must specify viscosity style for fix/rheo/viscosity"); +} + +/* ---------------------------------------------------------------------- */ + +FixRHEOViscosity::~FixRHEOViscosity() +{ + memory->destroy(eta_type); +} + +/* ---------------------------------------------------------------------- */ + +int FixRHEOViscosity::setmask() +{ + int mask = 0; + mask |= PRE_FORCE; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixRHEOViscosity::init() +{ + auto fixes = modify->get_fix_by_style("rheo"); + if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use fix rheo/viscosity"); + fix_rheo = dynamic_cast(fixes[0]); + + fix_rheo->viscosity_fix_defined = 1; + compute_grad = fix_rheo->compute_grad; +} + +/* ---------------------------------------------------------------------- */ + +void FixRHEOViscosity::setup_pre_force(int /*vflag*/) +{ + // Identify whether this is the last instance of fix viscosity + // Will handle communication + last_flag = 0; + + int i = 0; + auto fixlist = modify->get_fix_by_style("rheo/viscosity"); + for (const auto &ifix : fixlist) { + if (strcmp(ifix->id, id) == 0) break; + i++; + } + + if ((i + 1) == fixlist.size()) last_flag = 1; + + pre_force(0); +} + +/* ---------------------------------------------------------------------- */ + +void FixRHEOViscosity::pre_force(int /*vflag*/) +{ + int i, a, b; + double tmp, gdot; + + int *type = atom->type; + double *viscosity = fix_rheo->viscosity; + int *mask = atom->mask; + double **gradv = compute_grad->gradv; + + int nlocal = atom->nlocal; + int dim = domain->dimension; + + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + if (viscosity_style == CONSTANT) { + viscosity[i] = eta; + } else if (viscosity_style == TYPE) { + viscosity[i] = eta_type[type[i]]; + } else if (viscosity_style == POWER) { + gdot = 0.0; + for (a = 0; a < dim; a++) { + for (b = a; b < dim; b++) { + tmp = gradv[i][a * dim + b] + gradv[i][b * dim + a]; + tmp = tmp * tmp; + if (a == b) tmp *= 0.5; + gdot += tmp; + } + } + gdot = sqrt(gdot); + if (gdot <= gd0) { + viscosity[i] = eta; + } else { + viscosity[i] = K * pow(gdot, npow - 1) + tau0 / gdot; + } + } + } + } + + if (last_flag) comm->forward_comm(this); +} + +/* ---------------------------------------------------------------------- */ + +int FixRHEOViscosity::pack_forward_comm(int n, int *list, double *buf, + int /*pbc_flag*/, int * /*pbc*/) +{ + int i,j,k,m; + double *viscosity = fix_rheo->viscosity; + m = 0; + + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = viscosity[j]; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void FixRHEOViscosity::unpack_forward_comm(int n, int first, double *buf) +{ + int i, k, m, last; + double *viscosity = fix_rheo->viscosity; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + viscosity[i] = buf[m++]; + } +} diff --git a/src/RHEO/fix_rheo_viscosity.h b/src/RHEO/fix_rheo_viscosity.h new file mode 100644 index 0000000000..130d4d8cc8 --- /dev/null +++ b/src/RHEO/fix_rheo_viscosity.h @@ -0,0 +1,49 @@ +/* -*- 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/viscosity,FixRHEOViscosity) +// clang-format on +#else + +#ifndef LMP_FIX_RHEO_VISCOSITY_H +#define LMP_FIX_RHEO_VISCOSITY_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixRHEOViscosity : public Fix { + public: + FixRHEOViscosity(class LAMMPS *, int, char **); + ~FixRHEOViscosity(); + int setmask() override; + void init() override; + virtual void setup_pre_force(int) override; + void pre_force(int) override; + int pack_forward_comm(int, int *, double *, int, int *) override; + void unpack_forward_comm(int, int, double *) override; + private: + double *eta_type, eta; + double npow, K, gd0, tau0; + int viscosity_style; + int last_flag; + class FixRHEO *fix_rheo; + class ComputeRHEOGrad *compute_grad; +}; + +} // namespace LAMMPS_NS + +#endif +#endif