diff --git a/src/RHEO/compute_rheo_grad.cpp b/src/RHEO/compute_rheo_grad.cpp index f963485157..ce351048be 100644 --- a/src/RHEO/compute_rheo_grad.cpp +++ b/src/RHEO/compute_rheo_grad.cpp @@ -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; diff --git a/src/RHEO/compute_rheo_vshift.cpp b/src/RHEO/compute_rheo_vshift.cpp new file mode 100644 index 0000000000..42fa1e8c82 --- /dev/null +++ b/src/RHEO/compute_rheo_vshift.cpp @@ -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 +#include +#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; +} diff --git a/src/RHEO/compute_rheo_vshift.h b/src/RHEO/compute_rheo_vshift.h new file mode 100644 index 0000000000..fceb287510 --- /dev/null +++ b/src/RHEO/compute_rheo_vshift.h @@ -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 diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index d142af6483..cccaa792ff 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -132,7 +132,7 @@ FixRHEO::~FixRHEO() void FixRHEO::post_constructor() { - compute_kernel = dynamic_cast(modify->add_compute(fmt::format("rheo_kernel all rheo/kernel {} {} {}", kernel_style, zmin_kernel, h))); + compute_kernel = dynamic_cast(modify->add_compute("rheo_kernel all rheo/kernel")); fix_store_visc = dynamic_cast(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(modify->add_compute(fmt::format(cmd, h))); + compute_grad = dynamic_cast(modify->add_compute(cmd)); compute_grad->fix_rheo = this; if (rhosum_flag) - compute_rhosum = dynamic_cast(modify->add_compute(fmt::format("rheo_rhosum all rheo/rho/sum {} {}", h, zmin_rhosum))); + compute_rhosum = dynamic_cast(modify->add_compute("rheo_rhosum all rheo/rho/sum")); - if (shift_flag) - compute_vshift = dynamic_cast(modify->add_compute(fmt::format("rheo_vshift all rheo/vshift {}", h))); + if (shift_flag) { + compute_vshift = dynamic_cast(modify->add_compute("rheo_vshift all rheo/vshift")); + compute_vshift->fix_rheo = this; + } if (surface_flag) { fix_store_surf = dynamic_cast(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(modify->add_compute(fmt::format("rheo_interface all rheo/interface {}", h))); + compute_interface = dynamic_cast(modify->add_compute(fmt::format("rheo_interface all rheo/interface"))); fix_store_fp = dynamic_cast(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"); } /* ---------------------------------------------------------------------- */ diff --git a/src/RHEO/fix_rheo_pressure.cpp b/src/RHEO/fix_rheo_pressure.cpp new file mode 100644 index 0000000000..8f42e22239 --- /dev/null +++ b/src/RHEO/fix_rheo_pressure.cpp @@ -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 + +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(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; +} \ No newline at end of file diff --git a/src/RHEO/fix_rheo_pressure.h b/src/RHEO/fix_rheo_pressure.h new file mode 100644 index 0000000000..45217c7f75 --- /dev/null +++ b/src/RHEO/fix_rheo_pressure.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/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 diff --git a/src/RHEO/fix_rheo_viscosity.cpp b/src/RHEO/fix_rheo_viscosity.cpp index 3dfbfc1058..0fe9c7796e 100644 --- a/src/RHEO/fix_rheo_viscosity.cpp +++ b/src/RHEO/fix_rheo_viscosity.cpp @@ -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); diff --git a/src/RHEO/pair_rheo.cpp b/src/RHEO/pair_rheo.cpp index 06913e45e6..d3ac28d3d5 100644 --- a/src/RHEO/pair_rheo.cpp +++ b/src/RHEO/pair_rheo.cpp @@ -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(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(fixes[0]); + compute_kernel = fix_rheo->compute_kernel; compute_grad = fix_rheo->compute_grad; compute_interface = fix_rheo->compute_interface; diff --git a/src/RHEO/pair_rheo.h b/src/RHEO/pair_rheo.h index 18458048f6..aa98e63ddc 100644 --- a/src/RHEO/pair_rheo.h +++ b/src/RHEO/pair_rheo.h @@ -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