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

Drafting viscosity fix

Browse Source
This commit is contained in:
jtclemm
2023-02-19 22:08:02 -07:00
parent 6cbade576c
commit 560bd90e11
6 changed files with 272 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files

11
src/RHEO/compute_rheo_grad.cpp
View File

@ -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 <cmath>
#include <string>
@ -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++) {

16
src/RHEO/compute_rheo_grad.h
View File

@ -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;

9
src/RHEO/fix_rheo.cpp
View File

@ -123,6 +123,12 @@ void FixRHEO::post_constructor()
if (shift_flag)
compute_vshift = dynamic_cast<ComputeRHEOVShift *>(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");
}
/* ---------------------------------------------------------------------- */

3
src/RHEO/fix_rheo.h
View File

@ -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;

198
src/RHEO/fix_rheo_viscosity.cpp Normal file
View File

@ -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 <cmath>
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<FixRHEO *>(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++];
}
}

49
src/RHEO/fix_rheo_viscosity.h Normal file
View File

@ -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