Adding vshift + pressure, various fixes

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 <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;
+}

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

src/RHEO/fix_rheo.cpp

@@ -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)
+  if (shift_flag) {
-    compute_vshift = dynamic_cast<ComputeRHEOVShift *>(modify->add_compute(fmt::format("rheo_vshift all rheo/vshift {}", h)));
+    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");
 }
 
 /* ---------------------------------------------------------------------- */

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 <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;
+}

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

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;
@@ -69,33 +70,15 @@ void PairRHEO::compute(int eflag, int vflag)
   int error_flag, pair_force_flag, pair_rho_flag, pair_avisc_flag;
   double xtmp, ytmp, ztmp;
   int fluidi, fluidj;
-
+  double imass, jmass, rsq, r, rinv, psi_ij, Fij;
-  int *ilist, *jlist, *numneigh, **firstneigh;
-  double imass, jmass, rsq, r, rinv;
-
   double w, wp, rhoi, rhoj, voli, volj, Pi, Pj;
-  double *dWij, *dWji, *d2Wij, *d2Wji, *dW1ij, *dW1ji;
+  double *dWij, *dWji, *d2Wij, *d2Wji;
-  double vijeij, etai, etaj, kappai, kappaj;
+  double Ti, Tj, dT, etai, etaj, kappai, kappaj;
-  double Ti, Tj, dT;
+  double drho_damp, fmag, mu, q, fp_prefactor;
-  double drho_damp, fmag;
+  double dx[3], fv[3], dfp[3], fsolid[3], du[3], vi[3], vj[3];
-  double mu, q, fp_prefactor;
-  double dx[3] = {0};
-  double fv[3] = {0};
-  double dfp[3] = {0};
-  double fsolid[3] = {0};
-  double du[3] = {0};
-  double vi[3] = {0};
-  double vj[3] = {0};
-  double dv[3] = {0};
-  double psi_ij = 0.0;
-  double Fij = 0.0;
 
   ev_init(eflag, vflag);
 
-  double **gradv = compute_grad->gradv;
-  double **gradt = compute_grad->gradt;
-  double **gradr = compute_grad->gradr;
-  double **gradn = compute_grad->gradn;
   double **v = atom->v;
   double **x = atom->x;
   double **f = atom->f;
@@ -104,16 +87,20 @@ void PairRHEO::compute(int eflag, int vflag)
   double *drho = atom->drho;
   double *temperature = atom->temperature;
   double *heatflow = atom->heatflow;
+  double **gradv = compute_grad->gradv;
+  double **gradt = compute_grad->gradt;
+  double **gradr = compute_grad->gradr;
+  double **gradn = compute_grad->gradn;
   double *chi = compute_interface->chi;
   double **f_pressure = fix_rheo->f_pressure;
   double *viscosity = fix_rheo->viscosity;
   double *conductivity = fix_rheo->conductivity;
   double *pressure = fix_rheo->pressure;
-  double *special_lj = force->special_lj;
   tagint *tag = atom->tag;
   int *type = atom->type;
   int *status = atom->status;
 
+  int *ilist, *jlist, *numneigh, **firstneigh;
   int nlocal = atom->nlocal;
   int newton_pair = force->newton_pair;
   int dim = domain->dimension;
@@ -142,7 +129,6 @@ void PairRHEO::compute(int eflag, int vflag)
 
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
-      factor_lj = special_lj[sbmask(j)];
       j &= NEIGHMASK;
 
       dx[0] = xtmp - x[j][0];
@@ -152,7 +138,7 @@ void PairRHEO::compute(int eflag, int vflag)
       jtype = type[j];
       jmass = mass[jtype];
 
-      if (rsq < cutsq[itype][jtype]) {
+      if (rsq < hsq) {
         r = sqrt(rsq);
         rinv = 1 / r;
 
@@ -181,7 +167,6 @@ void PairRHEO::compute(int eflag, int vflag)
           vj[a] = v[j][a];
         }
 
-        // TODO: search for fix pressure and add calculate P function
         // Add corrections for walls
         rhoi = rho[i];
         rhoj = rho[j];
@@ -191,14 +176,14 @@ void PairRHEO::compute(int eflag, int vflag)
         if (fluidi && (!fluidj)) {
           compute_interface->correct_v(v[i], v[j], vi, i, j);
           rhoj = compute_interface->correct_rho(j, i);
-          Pj = fix_pressure->calculate_P(rhoj);
+          Pj = fix_rheo_pressure->calculate_p(rhoj);
           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(v[j], v[i], vj, j, i);
           rhoi = compute_interface->correct_rho(i,j);
-          Pi = fix_pressure->calculate_P(rhoi);
+          Pi = fix_rheo_pressure->calculate_p(rhoi);
           if ((chi[i] > 0.9) && (r < (h * 0.5)))
             fmag = (chi[i] - 0.9) * (h * 0.5 - r) * rho0 * csq * h * rinv;
 
@@ -429,9 +414,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;

src/RHEO/pair_rheo.h

@@ -50,6 +50,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