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

Fixing merge conflicts

Browse Source
This commit is contained in:
jtclemm
2023-03-22 21:35:54 -06:00
parent 886c642e01 6e65d13ad3
commit b550d72e38
9 changed files with 689 additions and 48 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/RHEO/compute_rheo_grad.cpp
View File

@ -39,7 +39,7 @@ ComputeRHEOGrad::ComputeRHEOGrad(LAMMPS *lmp, int narg, char **arg) :
if (narg < 4) error->all(FLERR,"Illegal compute rheo/grad command");
velocity_flag = temperature_flag = rho_flag = eta_flag = 0;
for (int iarg = 3; iarg < narg; iarg ++) {
for (int iarg = 3; iarg < narg; iarg++) {
if (strcmp(arg[iarg],"velocity") == 0) velocity_flag = 1;
else if (strcmp(arg[iarg],"rho") == 0) rho_flag = 1;
else if (strcmp(arg[iarg],"temperature") == 0) temperature_flag = 1;

294
src/RHEO/compute_rheo_vshift.cpp Normal file
View File

@ -0,0 +1,294 @@
/* ----------------------------------------------------------------------
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 "compute_rheo_vshift.h"
#include "fix_rheo.h"
#include "compute_rheo_solids.h"
#include "compute_rheo_grad.h"
#include "compute_rheo_kernel.h"
#include "fix_rheo_surface.h"
#include <cmath>
#include <cstring>
#include "atom.h"
#include "modify.h"
#include "domain.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "force.h"
#include "pair.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeRHEOVShift::ComputeRHEOVShift(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), vshift(nullptr), fix_rheo(nullptr), compute_kernel(nullptr),
compute_grad(nullptr), compute_surface(nullptr), compute_interface(nullptr)
{
if (narg != 3) error->all(FLERR,"Illegal compute RHEO/VShift command");
comm_reverse = 3;
surface_flag = 0;
nmax = atom->nmax;
memory->create(vshift, nmax, 3, "rheo/vshift:vshift");
array_atom = vshift;
peratom_flag = 1;
size_peratom_cols = 3;
}
/* ---------------------------------------------------------------------- */
ComputeRHEOVShift::~ComputeRHEOVShift()
{
memory->destroy(vshift);
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOVShift::init()
{
neighbor->add_request(this, NeighConst::REQ_DEFAULT);
surface_flag = 0;
if (fix_rheo->surface_flag) {
surface_flag = 1;
fix_rheo_surface = fix_rheo->fix_rheo_surface;
}
compute_kernel = fix_rheo->compute_kernel;
compute_grad = fix_rheo->compute_grad;
compute_interface = fix_rheo->compute_interface;
cut = fix_rheo->cut;
cutsq = cut * cut;
cutthird = cut / 3.0;
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOVShift::init_list(int /*id*/, NeighList *ptr)
{
list = ptr;
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOVShift::compute_peratom()
{
int i, j, a, b, ii, jj, jnum, itype, jtype;
int fluidi, fluidj;
double xtmp, ytmp, ztmp, rsq, r, rinv;
double w, wp, dr, w0, w4, vmag, prefactor;
double imass, jmass, voli, volj, rhoi, rhoj;
double dx[3], vi[3], vj[3] = {0};
int dim = domain->dimension;
int *jlist;
int inum, *ilist, *numneigh, **firstneigh;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
double **x = atom->x;
double **v = atom->v;
int *type = atom->type;
int *status = atom->status;
int *surface = atom->surface;
double *rho = atom->rho;
double *mass = atom->mass;
int newton_pair = force->newton_pair;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
if (nall > nmax) {
nmax = nall;
memory->destroy(vshift);
memory->create(vshift, nmax, 3, "rheo/vshift:vshift");
array_atom = vshift;
}
for (i = 0; i < nall; i++)
for (a = 0; a < dim; a++)
vshift[i][a] = 0.0;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
imass = mass[itype];
fluidi = status[i] & FixRHEO::STATUS_FLUID;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
fluidj = status[j] & FixRHEO::STATUS_FLUID;
if ((!fluidi) && (!fluidj)) continue;
if (!(status[i] & FixRHEO::STATUS_SHIFT) && !(status[j] & FixRHEO::STATUS_SHIFT)) continue;
dx[0] = xtmp - x[j][0];
dx[1] = ytmp - x[j][1];
dx[2] = ztmp - x[j][2];
rsq = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
if (rsq < cutsq) {
jtype = type[j];
jmass = mass[jtype];
r = sqrt(rsq);
rinv = 1 / r;
for (a = 0; a < dim; a ++) {
vi[a] = v[i][a];
vj[a] = v[j][a];
}
rhoi = rho[i];
rhoj = rho[j];
// Add corrections for walls
if (fluidi && (!fluidj)) {
compute_interface->correct_v(v[i], v[j], vi, i, j);
rhoj = compute_interface->correct_rho(j,i);
} else if ((!fluidi) && fluidj) {
compute_interface->correct_v(v[j], v[i], vj, j, i);
rhoi = compute_interface->correct_rho(i,j);
} else if ((!fluidi) && (!fluidj)) {
rhoi = 1.0;
rhoj = 1.0;
}
voli = imass / rhoi;
volj = jmass / rhoj;
wp = compute_kernel->calc_dw(i, j, dx[0], dx[1], dx[2], r);
w = compute_kernel->calc_w(i, j, dx[0], dx[1], dx[2], r);
w0 = compute_kernel->calc_w(i, j, 0, 0, 0, cutthird); // dx, dy, dz irrelevant
w4 = w * w * w * w / (w0 * w0 * w0 * w0);
dr = -2 * cutthird * (1 + 0.2 * w4) * wp * rinv;
if (mask[i] & groupbit) {
vmag = sqrt(vi[0] * vi[0] + vi[1] * vi[1] + vi[2] * vi[2]);
prefactor = vmag * volj * dr;
vshift[i][0] += prefactor * dx[0];
vshift[i][1] += prefactor * dx[1];
vshift[i][2] += prefactor * dx[2];
}
if (newton_pair || j < nlocal) {
if (mask[j] & groupbit) {
vmag = sqrt(vj[0] * vj[0] + vj[1] * vj[1] + vj[2] * vj[2]);
prefactor = vmag * voli * dr;
vshift[j][0] -= prefactor * dx[0];
vshift[j][1] -= prefactor * dx[1];
vshift[j][2] -= prefactor * dx[2];
}
}
}
}
}
if (newton_pair) comm->reverse_comm_compute(this);
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOVShift::correct_surfaces()
{
if (!surface_flag) return;
int *status = atom->status;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int i, a, b;
int dim = domain->dimension;
int *surface = atom->surface;
double **nsurf;
nsurf = fix_rheo_surface->n_surface;
double nx,ny,nz,vx,vy,vz;
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (surface[i] == 1 || surface[i] == 2) {
nx = nsurf[i][0];
ny = nsurf[i][1];
vx = vshift[i][0];
vy = vshift[i][1];
vz = vshift[i][2];
vshift[i][0] = (1 - nx * nx) * vx - nx * ny * vy;
vshift[i][1] = (1 - ny * ny) * vy - nx * ny * vx;
if (dim > 2) {
nz = nsurf[i][2];
vshift[i][0] -= nx * nz * vz;
vshift[i][1] -= ny * nz * vz;
vshift[i][2] = (1 - nz * nz) * vz - nz * ny * vy - nx * nz * vx;
} else {
vshift[i][2] = 0.0;
}
}
}
}
}
/* ---------------------------------------------------------------------- */
int ComputeRHEOVShift::pack_reverse_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
buf[m++] = vshift[i][0];
buf[m++] = vshift[i][1];
buf[m++] = vshift[i][2];
}
return m;
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOVShift::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
vshift[j][0] += buf[m++];
vshift[j][1] += buf[m++];
vshift[j][2] += buf[m++];
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeRHEOVShift::memory_usage()
{
double bytes = 3 * nmax * sizeof(double);
return bytes;
}

57
src/RHEO/compute_rheo_vshift.h Normal file
View File

@ -0,0 +1,57 @@
/* -*- 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/vshift,ComputeRHEOVShift)
// clang-format on
#else
#ifndef LMP_COMPUTE_RHEO_VSHIFT_H
#define LMP_COMPUTE_RHEO_VSHIFT_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeRHEOVShift : public Compute {
public:
ComputeRHEOVShift(class LAMMPS *, int, char **);
~ComputeRHEOVShift() override;
void init() override;
void init_list(int, class NeighList *) override;
void compute_peratom() override;
int pack_reverse_comm(int, int, double *) override;
void unpack_reverse_comm(int, int *, double *) override;
double memory_usage() override;
void correct_surfaces();
private:
int nmax;
double dtv, cut, cutsq, cutthird;
int surface_flag;
double **vshift;
class NeighList *list;
class FixRHEO *fix_rheo;
class FixRHEOSurface *fix_rheo_surface;
class ComputeRHEOInterface *compute_interface ;
class ComputeRHEOKernel *compute_kernel;
class ComputeRHEOGrad *compute_grad;
};
} // namespace LAMMPS_NS
#endif
#endif

45
src/RHEO/fix_rheo.cpp
View File

@ -132,7 +132,7 @@ FixRHEO::~FixRHEO()
void FixRHEO::post_constructor()
{
compute_kernel = dynamic_cast<ComputeRHEOKernel *>(modify->add_compute(fmt::format("rheo_kernel all rheo/kernel {} {} {}", kernel_style, zmin_kernel, h)));
compute_kernel = dynamic_cast<ComputeRHEOKernel *>(modify->add_compute("rheo_kernel all rheo/kernel"));
fix_store_visc = dynamic_cast<FixStorePeratom *>(modify->add_fix("rheo_store_visc all STORE/PERATOM 0 1"))
fix_store_visc->disable = 1;
@ -142,16 +142,18 @@ void FixRHEO::post_constructor()
fix_store_pres->disable = 1;
std::string cmd = "rheo_grad all rheo/grad {} velocity rho viscosity";
std::string cmd = "rheo_grad all rheo/grad velocity rho viscosity";
if (thermal_flag) cmd += "temperature";
compute_grad = dynamic_cast<ComputeRHEOGrad *>(modify->add_compute(fmt::format(cmd, h)));
compute_grad = dynamic_cast<ComputeRHEOGrad *>(modify->add_compute(cmd));
compute_grad->fix_rheo = this;
if (rhosum_flag)
compute_rhosum = dynamic_cast<ComputeRHEORhoSum *>(modify->add_compute(fmt::format("rheo_rhosum all rheo/rho/sum {} {}", h, zmin_rhosum)));
compute_rhosum = dynamic_cast<ComputeRHEORhoSum *>(modify->add_compute("rheo_rhosum all rheo/rho/sum"));
if (shift_flag)
compute_vshift = dynamic_cast<ComputeRHEOVShift *>(modify->add_compute(fmt::format("rheo_vshift all rheo/vshift {}", h)));
if (shift_flag) {
compute_vshift = dynamic_cast<ComputeRHEOVShift *>(modify->add_compute("rheo_vshift all rheo/vshift"));
compute_vshift->fix_rheo = this;
}
if (surface_flag) {
fix_store_surf = dynamic_cast<FixStorePeratom *>(modify->add_fix("rheo_store_surf all STORE/PERATOM 0 1"))
@ -166,7 +168,7 @@ void FixRHEO::post_constructor()
}
if (interface_flag) {
compute_interface = dynamic_cast<ComputeRHEOInterface *>(modify->add_compute(fmt::format("rheo_interface all rheo/interface {}", h)));
compute_interface = dynamic_cast<ComputeRHEOInterface *>(modify->add_compute(fmt::format("rheo_interface all rheo/interface")));
fix_store_fp = dynamic_cast<FixStorePeratom *>(modify->add_fix("rheo_store_fp all STORE/PERATOM 0 3"))
f_pressure = fix_store_fp->astore;
@ -212,7 +214,7 @@ void FixRHEO::setup_pre_force(int /*vflag*/)
void FixRHEO::setup()
{
// Confirm all accessory fixes are defined, may not cover group all
// Confirm all accessory fixes are defined
// Note: these fixes set this flag in setup_pre_force()
if (!viscosity_fix_defined)
error->all(FLERR, "Missing fix rheo/viscosity");
@ -231,6 +233,33 @@ void FixRHEO::setup()
viscosity_fix_defined = 0;
pressure_fix_defined = 0;
surface_fix_defined = 0;
// Check fixes cover all atoms (doesnt ensure user covers atoms created midrun)
// (pressure is currently required to be group all)
auto visc_fixes = modify->get_fix_by_style("rheo/viscosity");
auto therm_fixes = modify->get_fix_by_style("rheo/thermal");
int *mask = atom->mask;
int v_coverage_flag = 1;
int t_coverage_flag = 1;
int covered;
for (int i = 0; i < atom->nlocal; i++) {
covered = 0;
for (auto fix in visc_fixes)
if (mask[i] & fix->groupbit) covered = 1;
if (!covered) v_coverage_flag = 0;
if (thermal_flag) {
covered = 0;
for (auto fix in therm_fixes)
if (mask[i] & fix->groupbit) covered = 1;
if (!covered) v_coverage_flag = 0;
}
}
if (!v_coverage_flag)
error->one(FLERR, "Fix rheo/viscosity does not fully cover all atoms");
if (!t_coverage_flag)
error->one(FLERR, "Fix rheo/thermal does not fully cover all atoms");
}
/* ---------------------------------------------------------------------- */

206
src/RHEO/fix_rheo_pressure.cpp Normal file
View File

@ -0,0 +1,206 @@
/* ----------------------------------------------------------------------
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_pressure.h"
#include "atom.h"
#include "comm.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, LINEAR, CUBIC, TAITWATER};
/* ---------------------------------------------------------------------- */
FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 4) error->all(FLERR,"Illegal fix command");
pressure_style = NONE;
comm_forward = 1;
nmax = atom->nmax;
// Currently can only have one instance of fix rheo/pressure
if (igroup != 0)
error->all(FLERR,"fix rheo/pressure command requires group all");
int ntypes = atom->ntypes;
int iarg = 3;
if (strcmp(arg[iarg],"linear") == 0) {
pressure_style = LINEAR;
} else if (strcmp(arg[iarg],"taitwater") == 0) {
pressure_style = TAITWATER;
} else if (strcmp(arg[iarg],"cubic") == 0) {
pressure_style = CUBIC;
if (iarg + 1 >= narg) error->all(FLERR,"Insufficient arguments for pressure option");
c_cubic = utils::numeric(FLERR,arg[iarg + 1],false,lmp);
} else {
error->all(FLERR,"Illegal fix command, {}", arg[iarg]);
}
if (pressure_style == NONE)
error->all(FLERR,"Must specify pressure style for fix/rheo/pressure");
}
/* ---------------------------------------------------------------------- */
FixRHEOPressure::~FixRHEOPressure() {}
/* ---------------------------------------------------------------------- */
int FixRHEOPressure::setmask()
{
int mask = 0;
mask |= PRE_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixRHEOPressure::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/pressure");
fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
csq = fix_rheo->csq;
rho0 = fix_rheo->rho0;
rho0inv = 1.0 / rho0;
// Cannot define multiple as pair rheo cannot currently distinguish
if (modify->get_fix_by_style("rheo/pressure").size() > 1)
error->all(FLERR,"Can only specify one instance of fix rheo/pressure");
}
/* ---------------------------------------------------------------------- */
void FixRHEOPressure::setup_pre_force(int /*vflag*/)
{
fix_rheo->pressure_fix_defined = 1;
// Identify whether this is the first/last instance of fix pressure
// First will handle growing arrays
// Last will handle communication
first_flag = 0
last_flag = 0;
int i = 0;
auto fixlist = modify->get_fix_by_style("rheo/pressure");
for (const auto &ifix : fixlist) {
if (strcmp(ifix->id, id) == 0) break;
i++;
}
if (i == 0) first_flag = 1;
if ((i + 1) == fixlist.size()) last_flag = 1;
pre_force(0);
}
/* ----------------------------------------------------------------------
Update (and forward) pressure every timestep
------------------------------------------------------------------------- */
void FixRHEOPressure::pre_force(int /*vflag*/)
{
int i;
double dr, rr3, rho_ratio;
double *p = fix_rheo->pressure;
int *mask = atom->mask;
double *rho = atom->rho;
int nlocal = atom->nlocal;
int dim = domain->dimension;
if (first_flag & nmax < atom->nmax) {
nmax = atom->nmax;
fix_rheo->fix_store_visc->grow_arrays(nmax);
}
if (pressure_style == TAITWATER) inv7 = 1.0 / 7.0;
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (pressure_style == LINEAR) {
p[i] = csq * (rho[i] - rho0);
} else if (pressure_style == CUBIC) {
dr = rho[i] - rho0;
p[i] = csq * (dr + c_cubic * dr * dr * dr);
} else if (pressure_style == TAITWATER) {
rho_ratio = rho[i] / rho0inv;
rr3 = rho_ratio * rho_ratio * rho_ratio;
p[i] = csq * rho0 * inv7 * (rr3 * rr3 * rho_ratio - 1.0);
}
}
}
if (last_flag && comm_forward) comm->forward_comm(this);
}
/* ---------------------------------------------------------------------- */
int FixRHEOPressure::pack_forward_comm(int n, int *list, double *buf,
int /*pbc_flag*/, int * /*pbc*/)
{
int i,j,k,m;
double *pressure = fix_rheo->pressure;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = pressure[j];
}
return m;
}
/* ---------------------------------------------------------------------- */
void FixRHEOPressure::unpack_forward_comm(int n, int first, double *buf)
{
int i, k, m, last;
double *pressure = fix_rheo->pressure;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
pressure[i] = buf[m++];
}
}
double FixRHEOPressure::calculate_p(double rho)
{
double rho;
if (pressure_style == LINEAR) {
p = csq * (rho - rho0);
} else if (pressure_style == CUBIC) {
dr = rho - rho0;
p = csq * (dr + c_cubic * dr * dr * dr);
} else if (pressure_style == TAITWATER) {
rho_ratio = rho / rho0inv;
rr3 = rho_ratio * rho_ratio * rho_ratio;
p = csq * rho0 * inv7 * (rr3 * rr3 * rho_ratio - 1.0);
}
return rho;
}

49
src/RHEO/fix_rheo_pressure.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/pressure,FixRHEOPressure)
// clang-format on
#else
#ifndef LMP_FIX_RHEO_PRESSURE_H
#define LMP_FIX_RHEO_PRESSURE_H
#include "fix.h"
namespace LAMMPS_NS {
class FixRHEOPressure : public Fix {
public:
FixRHEOPressure(class LAMMPS *, int, char **);
~FixRHEOPressure() override;
int setmask() override;
void init() override;
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;
double calculate_p(double);
private:
double c_cubic, csq, rho0, rho0inv;
int pressure_style;
int first_flag, last_flag;
int nmax;
class FixRHEO *fix_rheo;
};
} // namespace LAMMPS_NS
#endif
#endif

6
src/RHEO/fix_rheo_viscosity.cpp
View File

@ -44,13 +44,13 @@ FixRHEOViscosity::FixRHEOViscosity(LAMMPS *lmp, int narg, char **arg) :
int ntypes = atom->ntypes;
int iarg = 3;
if (strcmp(arg[iarg],"constant") == 0) {
if (iarg+1 >= narg) error->all(FLERR,"Insufficient arguments for viscosity option");
if (iarg + 1 >= narg) error->all(FLERR,"Insufficient arguments for viscosity option");
viscosity_style = CONSTANT;
eta = utils::numeric(FLERR,arg[iarg + 1],false,lmp);
if (eta < 0.0) error->all(FLERR,"The viscosity must be positive");
iarg += 1;
} else if (strcmp(arg[iarg],"type") == 0) {
if(iarg+ntypes >= narg) error->all(FLERR,"Insufficient arguments for viscosity option");
if (iarg + ntypes >= narg) error->all(FLERR,"Insufficient arguments for viscosity option");
viscosity_style = TYPE;
memory->create(eta_type, ntypes + 1, "rheo_thermal:eta_type");
for (int i = 1; i <= ntypes; i++) {
@ -59,7 +59,7 @@ FixRHEOViscosity::FixRHEOViscosity(LAMMPS *lmp, int narg, char **arg) :
}
iarg += ntypes;
} else if (strcmp(arg[iarg],"power") == 0) {
if (iarg+4 >= narg) error->all(FLERR,"Insufficient arguments for viscosity option");
if (iarg + 4 >= narg) error->all(FLERR,"Insufficient arguments for viscosity option");
viscosity_style = POWER;
comm_forward = 1;
eta = utils::numeric(FLERR,arg[iarg + 1],false,lmp);

76
src/RHEO/pair_rheo.cpp
View File

@ -21,6 +21,7 @@
#include "domain.h"
#include "error.h"
#include "fix_rheo.h"
#include "fix_rheo_pressure.h"
#include "force.h"
#include "math_extra.h"
#include "memory.h"
@ -41,7 +42,7 @@ using namespace MathExtra;
PairRHEO::PairRHEO(LAMMPS *lmp) :
Pair(lmp), compute_kernel(nullptr), compute_grad(nullptr),
compute_interface(nullptr), fix_rheo(nullptr)
compute_interface(nullptr), fix_rheo(nullptr), fix_rheo_pressure(nullptr)
{
restartinfo = 0;
single_enable = 0;
@ -67,6 +68,7 @@ void PairRHEO::compute(int eflag, int vflag)
{
int i, j, a, b, ii, jj, inum, jnum, itype, jtype;
int pair_force_flag, pair_rho_flag, pair_avisc_flag;
int fluidi, fluidj;
double xtmp, ytmp, ztmp, w, wp, Ti, Tj, dT;
double rhoi, rhoj, voli, volj, Pi, Pj, etai, etaj, kappai, kappaj;
double mu, q, fp_prefactor, drho_damp, fmag, psi_ij, Fij;
@ -84,15 +86,16 @@ void PairRHEO::compute(int eflag, int vflag)
double **v = atom->v;
double **x = atom->x;
double **f = atom->f;
double **fp = atom->fp; // rewrite later
double **f_pressure = fix_rheo->f_pressure; // rewrite later
double *pressure = atom->pressure; // rewrite later
double *rho = atom->rho;
double *mass = atom->mass;
double *drho = atom->drho;
double *temp = atom->temp;
double *heat = atom->heat;
double *temperature = atom->temperature;
double *heatflow = atom->heatflow;
double *special_lj = force->special_lj;
tagint *tag = atom->tag;
int *chi = compute_interface->chi;
int *type = atom->type;
int *status = atom->status;
@ -108,6 +111,7 @@ void PairRHEO::compute(int eflag, int vflag)
conductivity = atom->dvector[index_cond];
}
int *ilist, *jlist, *numneigh, **firstneigh;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
int dim = domain->dimension;
@ -129,14 +133,14 @@ void PairRHEO::compute(int eflag, int vflag)
jnum = numneigh[i];
imass = mass[itype];
etai = viscosity[i];
fluidi = status[i] & FixRHEO::STATUS_FLUID;
if (thermal_flag) {
kappai = conductivity[i];
Ti = temp[i];
Ti = temperature[i];
}
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
j &= NEIGHMASK;
dx[0] = xtmp - x[j][0];
@ -145,24 +149,25 @@ void PairRHEO::compute(int eflag, int vflag)
rsq = lensq3(dx);
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
if (rsq < hsq) {
r = sqrt(rsq);
rinv = 1 / r;
jmass = mass[jtype];
etaj = viscosity[j];
fluidj = status[j] & FixRHEO::STATUS_FLUID;
if (thermal_flag) {
Tj = temp[j];
Tj = temperature[j];
kappaj = conductivity[j];
}
pair_rho_flag = 0;
pair_force_flag = 0;
pair_avisc_flag = 0;
if (status[i] <= FixRHEO::FLUID_MAX || status[j] <= FixRHEO::FLUID_MAX) {
if (fluidi || fluidj) {
pair_force_flag = 1;
}
if (status[i] <= FixRHEO::FLUID_MAX && status[j] <= FixRHEO::FLUID_MAX) {
if (fluidi && fluidj) {
pair_avisc_flag = 1;
pair_rho_flag = 1;
}
@ -171,10 +176,9 @@ void PairRHEO::compute(int eflag, int vflag)
dWij = compute_kernel->dWij;
dWji = compute_kernel->dWji;
for (a = 0; a < dim; a ++) {
for (a = 0; a < dim; a++) {
vi[a] = v[i][a];
vj[a] = v[j][a];
fsolid[a] = 0.0;
}
// Add corrections for walls
@ -182,48 +186,46 @@ void PairRHEO::compute(int eflag, int vflag)
rhoj = rho[j];
Pi = pressure[i];
Pj = pressure[j];
if ((status[i] & STATUS_FLUID) && !(status[j] & STATUS_FLUID)) {
fmag = 0;
if (fluidi && (!fluidj)) {
compute_interface->correct_v(vi, vj, i, j);
rhoj = compute_interface->correct_rho(j, i);
Pj = calc_pressure(rhoj, jtype);
Pj = fix_rheo_pressure->calculate_p(rhoj);
// Repel if close to inner solid particle
if (compute_interface->chi[j] > 0.9 && r < (h * 0.5)) {
fmag = (compute_interface->chi[j] - 0.9) * (h * 0.5 - r);
fmag *= rho0 * csq * h * ir;
scale3(fmag, dx, fsolid);
}
} else if (!(status[i] & STATUS_FLUID) && (status[j] & STATUS_FLUID)) {
if ((chi[j] > 0.9) && (r < (h * 0.5)))
fmag = (chi[j] - 0.9) * (h * 0.5 - r) * rho0 * csq * h * rinv;
} else if ((!fluidi) && fluidj) {
compute_interface->correct_v(vj, vi, j, i);
rhoi = compute_interface->correct_rho(i, j);
Pi = calc_pressure(rhoi, itype);
// Repel if close to inner solid particle
if (compute_interface->chi[i] > 0.9 && r < (h * 0.5)) {
fmag = (compute_interface->chi[i] - 0.9) * (h * 0.5 - r);
fmag *= rho0 * csq * h * ir;
scale3(fmag, dx, fsolid);
}
} else if (!(status[i] & STATUS_FLUID) && !(status[j] & STATUS_FLUID)) {
if (chi[i] > 0.9 && r < (h * 0.5)) {
fmag = (chi[i] - 0.9) * (h * 0.5 - r) * rho0 * csq * h * rinv;
} else if ((!fluidi) && (!fluidj)) {
rhoi = 1.0;
rhoj = 1.0;
}
// Compute volume and pressure after reconstructing
// Repel if close to inner solid particle
scale3(fmag, dx, fsolid);
// Compute volume after reconstructing
voli = imass / rhoi;
volj = jmass / rhoj;
//Thermal Evolution
// Thermal Evolution
if (thermal_flag) {
dT = dot3(dx, dWij);
dT *= (kappai + kappaj) * (Ti - Tj) * rinv * rinv * voli * volj;
//Assumes heat capacity and density = 1, needs to be generalized
heat[i] += dT;
//TODO: Assumes heat capacity and density = 1, needs to be generalized
heatflow[i] += dT;
if (newton_pair || j < nlocal) {
dT = dot3(dx, dWji);
dT *= (kappai + kappaj) * (Tj - Ti) * rinv * rinv * voli * volj;
heat[j] -= dT;
heatflow[j] -= dT;
}
}
@ -242,7 +244,7 @@ void PairRHEO::compute(int eflag, int vflag)
du[a] -= 0.5 * (gradv[i][a * dim + b] + gradv[j][a * dim + b]) * dx[b];
mu = dot3(du, dx) * hinv3;
mu = mu / (rsq * hinv3 * hinv3 + EPSILON);
mu /= (rsq * hinv3 * hinv3 + EPSILON);
mu = MIN(0.0, mu);
q = av * (-2.0 * cs * mu + mu * mu);
fp_prefactor += voli * volj * q * (rhoj + rhoi);
@ -403,9 +405,13 @@ void PairRHEO::coeff(int narg, char **arg)
void PairRHEO::setup()
{
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");
if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use pair rheo");
fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
fixes = modify->get_fix_by_style("rheo/pressure");
if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo/pressure to use pair rheo");
fix_rheo_pressure = dynamic_cast<FixRHEOPressure *>(fixes[0]);
compute_kernel = fix_rheo->compute_kernel;
compute_grad = fix_rheo->compute_grad;
compute_interface = fix_rheo->compute_interface;

2
src/RHEO/pair_rheo.h
View File

@ -36,7 +36,6 @@ class PairRHEO : public Pair {
protected:
double h, csq, rho0; // From fix RHEO
double cs, hsq, hinv, hinv3, av, rho_damp;
int laplacian_order;
@ -50,6 +49,7 @@ class PairRHEO : public Pair {
class ComputeRHEOGrad *compute_grad;
class ComputeRHEOInterface *compute_interface;
class FixRHEO *fix_rheo;
class FixRHEOPressure *fix_rheo_pressure;
};
} // namespace LAMMPS_NS