merge conflicts

src/.gitignore

@@ -217,6 +217,8 @@
 /fix_rheo.h
 /fix_rheo_pressure.cpp
 /fix_rheo_pressure.h
+/fix_rheo_stress.cpp
+/fix_rheo_stress.h
 /fix_rheo_thermal.cpp
 /fix_rheo_thermal.h
 /fix_rheo_viscosity.cpp
@@ -225,6 +227,8 @@
 /pair_rheo.h
 /pair_rheo_react.cpp
 /pair_rheo_react.h
+/pair_rheo_tension.cpp
+/pair_rheo_tension.h
 
 /compute_grid.cpp
 /compute_grid.h

src/RHEO/compute_rheo_grad.cpp

@@ -111,6 +111,8 @@ void ComputeRHEOGrad::init()
   compute_kernel = fix_rheo->compute_kernel;
   compute_interface = fix_rheo->compute_interface;
 
+  remap_v_flag = domain->deform_vremap;
+
   neighbor->add_request(this, NeighConst::REQ_DEFAULT);
 }
 
@@ -301,13 +303,23 @@ void ComputeRHEOGrad::forward_fields()
 /* ---------------------------------------------------------------------- */
 
 int ComputeRHEOGrad::pack_forward_comm(int n, int *list, double *buf,
-                                        int /*pbc_flag*/, int * /*pbc*/)
+                                        int pbc_flag, int *pbc)
 {
   int i,j,k,m;
+  int *mask = atom->mask;
   double *rho = atom->rho;
   double *temperature = atom->temperature;
   double **v = atom->v;
   int dim = domain->dimension;
+  double *h_rate = domain->h_rate;
+  int deform_groupbit = domain->deform_groupbit;
+  double dv[3];
+
+  if (remap_v_flag) {
+    dv[0] = pbc[0] * h_rate[0] + pbc[5] * h_rate[5] + pbc[4] * h_rate[4];
+    dv[1] = pbc[1] * h_rate[1] + pbc[3] * h_rate[3];
+    dv[2] = pbc[2] * h_rate[2];
+  }
 
   m = 0;
 
@@ -333,9 +345,17 @@ int ComputeRHEOGrad::pack_forward_comm(int n, int *list, double *buf,
 
     } else if (comm_stage == COMMFIELD) {
 
-      if (velocity_flag)
-        for (k = 0; k < dim; k++)
-          buf[m++] = v[j][k];
+      if (velocity_flag) {
+        if (remap_v_flag & pbc_flag & (mask[j] & deform_groupbit)) {
+          for (k = 0; k < dim; k++)
+            buf[m++] = v[j][k] + dv[k];
+        } else {
+          for (k = 0; k < dim; k++)
+            buf[m++] = v[j][k];
+        }
+      }
+
+
 
       if (rho_flag)
         buf[m++] = rho[j];

src/RHEO/compute_rheo_grad.h

@@ -49,7 +49,7 @@ class ComputeRHEOGrad : public Compute {
   double cut, cutsq, rho0;
 
   int velocity_flag, temperature_flag, rho_flag, eta_flag;
-  int interface_flag;
+  int interface_flag, remap_v_flag;
 
   class ComputeRHEOKernel *compute_kernel;
   class ComputeRHEOInterface *compute_interface;

src/RHEO/fix_rheo.cpp

@@ -381,13 +381,13 @@ void FixRHEO::pre_force(int /*vflag*/)
   if (rhosum_flag)
     compute_rhosum->compute_peratom();
 
-  compute_grad->forward_fields(); // also forwards v and rho for chi
   compute_kernel->compute_peratom();
   if (interface_flag) {
     // Note on first setup, have no forces for pressure to reference
     compute_interface->compute_peratom();
   }
 
+  // No need to forward v, rho, or T for compute_grad since already done
   compute_grad->compute_peratom();
   compute_grad->forward_gradients();
 

src/RHEO/fix_rheo_stress.cpp

@@ -0,0 +1,137 @@
+/* ----------------------------------------------------------------------
+ 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.
+ ------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing authors: Joel Clemmer (SNL)
+----------------------------------------------------------------------- */
+
+#include "fix_rheo_stress.h"
+
+#include "atom.h"
+#include "comm.h"
+#include "compute.h"
+#include "domain.h"
+#include "fix_store_atom.h"
+#include "group.h"
+#include "error.h"
+#include "modify.h"
+#include "update.h"
+
+#include <cmath>
+#include <cstring>
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+/* ---------------------------------------------------------------------- */
+
+FixRHEOStress::FixRHEOStress(LAMMPS *lmp, int narg, char **arg) :
+  id_compute(nullptr), id_fix(nullptr), stress_compute(nullptr), store_fix(nullptr), Fix(lmp, narg, arg)
+{
+  if (narg != 3) error->all(FLERR,"Illegal fix rheo/stress command");
+  comm_forward = 6;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixRHEOStress::~FixRHEOStress()
+{
+  modify->delete_compute(id_compute);
+  modify->delete_fix(id_fix);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixRHEOStress::post_constructor()
+{
+  id_fix = utils::strdup(std::string(id) + "_store");
+  store_fix = dynamic_cast<FixStoreAtom *>(modify->add_fix(fmt::format("{} {} STORE/ATOM d_pxx d_pyy d_pzz d_pxy d_pxz d_pyz", id_fix, group->names[igroup])));
+  array_atom = store_fix->astore;
+
+  id_compute = utils::strdup(std::string(id) + "_compute");
+  stress_compute = modify->add_compute(fmt::format("{} {} stress/atom NULL ke pair bond", id_compute, group->names[igroup]));
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixRHEOStress::setmask()
+{
+  int mask = 0;
+  mask |= PRE_FORCE;
+  mask |= END_OF_STEP;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixRHEOStress::init()
+{
+  stress_compute->addstep(update->ntimestep+1);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixRHEOStress::pre_force(int vflag)
+{
+  // add pre-force and forward to ghosts (not done in store/atom)
+  comm->forward_comm(this);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixRHEOStress::end_of_step()
+{
+  stress_compute->compute_peratom();
+
+  // copy compute to fix property atom
+  double **saved_stress = store_fix->astore;
+  double **stress = stress_compute->array_atom;
+
+  int ntotal = atom->nlocal+atom->nghost;
+  for (int i = 0; i < ntotal; i++)
+    for (int a = 0; a < 6; a++)
+      saved_stress[i][a] = stress[i][a];
+
+  stress_compute->addstep(update->ntimestep + 1);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixRHEOStress::pack_forward_comm(int n, int *list, double *buf,
+                                        int /*pbc_flag*/, int * /*pbc*/)
+{
+  int i, j, a, m;
+  double **saved_stress = store_fix->astore;
+
+  m = 0;
+  for (i = 0; i < n; i++) {
+    j = list[i];
+    for (a = 0; a < 6; a++)
+      buf[m++] = saved_stress[j][a];
+  }
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixRHEOStress::unpack_forward_comm(int n, int first, double *buf)
+{
+  int i, a, m, last;
+  double **saved_stress = store_fix->astore;
+
+  m = 0;
+  last = first + n;
+  for (i = first; i < last; i++)
+    for (a = 0; a < 6; a++)
+      saved_stress[i][a] = buf[m++];
+}

src/RHEO/fix_rheo_stress.h

@@ -0,0 +1,48 @@
+/* -*- 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/stress,FixRHEOStress);
+// clang-format on
+#else
+
+#ifndef LMP_FIX_RHEO_STRESS_H
+#define LMP_FIX_RHEO_STRESS_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixRHEOStress : public Fix {
+ public:
+  FixRHEOStress(class LAMMPS *, int, char **);
+  ~FixRHEOStress() override;
+  void post_constructor() override;
+  int setmask() override;
+  void init() override;
+  void pre_force(int) override;
+  void end_of_step() override;
+  int pack_forward_comm(int, int *, double *, int, int *) override;
+  void unpack_forward_comm(int, int, double *) override;
+
+ private:
+  char *id_compute, *id_fix;
+  class Compute *stress_compute;
+  class FixStoreAtom *store_fix;
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif

src/RHEO/pair_rheo_tension.cpp

@@ -0,0 +1,365 @@
+/* ----------------------------------------------------------------------
+ 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.
+ ------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing authors:
+   Joel Clemmer (SNL)
+----------------------------------------------------------------------- */
+
+#include "pair_rheo_tension.h"
+
+#include "atom.h"
+#include "comm.h"
+#include "compute_rheo_kernel.h"
+#include "domain.h"
+#include "error.h"
+#include "fix_rheo.h"
+#include "force.h"
+#include "math_extra.h"
+#include "memory.h"
+#include "modify.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "update.h"
+#include "utils.h"
+
+#include <cmath>
+
+using namespace LAMMPS_NS;
+using namespace RHEO_NS;
+using namespace MathExtra;
+
+static constexpr double EPSILON = 1e-2;
+
+/* ---------------------------------------------------------------------- */
+
+PairRHEOTension::PairRHEOTension(LAMMPS *lmp) :
+  Pair(lmp), compute_kernel(nullptr), fix_rheo(nullptr)
+{
+  restartinfo = 0;
+  single_enable = 0;
+
+  comm_forward = 3;
+  comm_reverse = 3;
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairRHEOTension::~PairRHEOTension()
+{
+  // Remove custom property if it exists
+  int tmp1, tmp2, index;
+
+  index = atom->find_custom("rheo_c_tension", tmp1, tmp2);
+  if (index != -1) atom->remove_custom(index, 1, 0);
+
+  index = atom->find_custom("rheo_n_tension", tmp1, tmp2);
+  if (index != -1) atom->remove_custom(index, 1, 3);
+
+  if (allocated) {
+    memory->destroy(alpha);
+    memory->destroy(setflag);
+    memory->destroy(cutsq);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairRHEOTension::compute(int eflag, int vflag)
+{
+  int i, j, a, b, ii, jj, inum, jnum, itype, jtype;
+  int fluidi, fluidj;
+  double xtmp, ytmp, ztmp, w, wp;
+  double rhoi, rhoj, voli, volj;
+  double *dWij, *dWji;
+  double dx[3], ft[3];
+
+  int *ilist, *jlist, *numneigh, **firstneigh;
+  double imass, jmass, rsq, r, rinv;
+
+  int nlocal = atom->nlocal;
+  int newton_pair = force->newton_pair;
+  int dim = domain->dimension;
+
+  ev_init(eflag, vflag);
+
+  double **x = atom->x;
+  double **f = atom->f;
+  double *rho = atom->rho;
+  double *mass = atom->mass;
+  double *special_lj = force->special_lj;
+  int *type = atom->type;
+  int *status = atom->status;
+  tagint *tag = atom->tag;
+
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  int nmax = atom->nmax;
+  if (nmax_store <= nmax) {
+    memory->grow(ct, nmax, "atom:rheo_c_tension");
+    memory->grow(nnt_tension, nmax, 3, "atom:rheo_n_tension");
+    nmax_store = atom->nmax;
+  }
+
+  // loop over neighbors of my atoms
+
+  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];
+    rhoi = rho[i];
+    voli = imass / rhoi;
+    fluidi = !(status[i] & PHASECHECK);
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      j &= NEIGHMASK;
+
+      dx[0] = xtmp - x[j][0];
+      dx[1] = ytmp - x[j][1];
+      dx[2] = ztmp - x[j][2];
+      rsq = lensq3(dx);
+      jtype = type[j];
+
+      if (rsq > hsq) continue;
+
+      r = sqrt(rsq);
+      rinv = 1 / r;
+
+      jmass = mass[jtype];
+      rhoj = rho[j];
+      volj = jmass / rhoj;
+      fluidj = !(status[j] & PHASECHECK);
+
+      wp = compute_kernel->calc_dw(i, j, dx[0], dx[1], dx[2],r);
+      dWij = compute_kernel->dWij;
+      dWji = compute_kernel->dWji;
+
+      f[i][0] += ft[0];
+      f[i][1] += ft[1];
+      f[i][2] += ft[2];
+
+      if (evflag) // Does not account for unbalanced forces
+        ev_tally_xyz(i, j, nlocal, newton_pair, 0.0, 0.0, ft[0], ft[1], ft[2], dx[0], dx[1], dx[2]);
+
+      if (newton_pair || j < nlocal) {
+
+        f[j][0] -= ft[0];
+        f[j][1] -= ft[1];
+        f[j][2] -= ft[2];
+      }
+    }
+  }
+
+  if (vflag_fdotr) virial_fdotr_compute();
+
+  comm->reverse_comm(this);
+  comm->forward_comm(this);
+}
+
+/* ----------------------------------------------------------------------
+ allocate all arrays
+ ------------------------------------------------------------------------- */
+
+void PairRHEOTension::allocate()
+{
+  allocated = 1;
+  int n = atom->ntypes;
+
+  memory->create(setflag, n + 1, n + 1, "pair:setflag");
+  for (int i = 1; i <= n; i++)
+    for (int j = i; j <= n; j++)
+      setflag[i][j] = 0;
+
+  memory->create(alpha, n + 1, n + 1, "pair:alpha");
+  memory->create(cutsq, n + 1, n + 1, "pair:cutsq");
+}
+
+/* ----------------------------------------------------------------------
+ set coeffs for one or more type pairs
+ ------------------------------------------------------------------------- */
+
+void PairRHEOTension::coeff(int narg, char **arg)
+{
+  if (narg != 3)
+    error->all(FLERR,"Incorrect number of args for pair_style rheo coefficients");
+  if (!allocated)
+    allocate();
+
+  int ilo, ihi, jlo, jhi;
+  utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi,error);
+  utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi,error);
+
+  alpha_one = utils::numeric(FLERR, arg[2], false, lmp);
+
+  int count = 0;
+  for (int i = ilo; i <= ihi; i++) {
+    for (int j = 0; j <= atom->ntypes; j++) {
+      alpha[i][j] = alpha_one;
+      setflag[i][j] = 1;
+      count++;
+    }
+  }
+
+  if (count == 0)
+    error->all(FLERR,"Incorrect args for pair rheo/tension coefficients");
+}
+
+/* ----------------------------------------------------------------------
+ setup specific to this pair style
+ ------------------------------------------------------------------------- */
+
+void PairRHEOTension::setup()
+{
+  auto fixes = modify->get_fix_by_style("rheo");
+  if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use pair rheo");
+  fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
+
+  compute_kernel = fix_rheo->compute_kernel;
+  compute_grad = fix_rheo->compute_grad;
+  compute_interface = fix_rheo->compute_interface;
+  h = fix_rheo->h;
+  csq = fix_rheo->csq;
+  rho0 = fix_rheo->rho0;
+
+  hsq = h * h;
+  hinv = 1.0 / h;
+  hinv3 = hinv * 3.0;
+}
+
+/* ----------------------------------------------------------------------
+   init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairRHEOTension::init_style()
+{
+  neighbor->add_request(this);
+
+
+  // Create c_tension arrays n_tension arrays if they don't already exist
+  // Create a custom atom property so it works with compute property/atom
+  // Do not create grow callback as there's no reason to copy/exchange data
+  // Manually grow if nmax_store exceeded
+  // For B and gradC, create a local array since they are unlikely to be printed
+
+  int tmp1, tmp2;
+  int index = atom->find_custom("rheo_c_tension", tmp1, tmp2);
+  if (index == -1)  index = atom->add_custom("rheo_c_tension", 1, 0);
+  ct = atom->dvector[index];
+
+  index = atom->find_custom("rheo_n_tension", tmp1, tmp2);
+  if (index == -1)  index = atom->add_custom("rheo_n_tension", 1, 3);
+  nt = atom->darray[index];
+
+  nmax_store = atom->nmax;
+}
+
+/* ----------------------------------------------------------------------
+ init for one type pair i,j and corresponding j,i
+ ------------------------------------------------------------------------- */
+
+double PairRHEOTension::init_one(int i, int j)
+{
+  if (setflag[i][j] == 0)
+    error->all(FLERR,"All pair rheo/tension coeffs are not set");
+
+  alpha[j][i] = alpha[i][j];
+
+  return h;
+}
+
+
+/* ---------------------------------------------------------------------- */
+
+int PairRHEOTension::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/)
+{
+  int i,j,k,m;
+  m = 0;
+  double *rho = atom->rho;
+
+  for (i = 0; i < n; i++) {
+    j = list[i];
+    if (comm_stage == 0) {
+      buf[m++] = fp_store[j][0];
+      buf[m++] = fp_store[j][1];
+      buf[m++] = fp_store[j][2];
+    } else {
+      buf[m++] = chi[j];
+      buf[m++] = rho[j];
+    }
+  }
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairRHEOTension::unpack_forward_comm(int n, int first, double *buf)
+{
+  int i, k, m, last;
+  double *rho = atom->rho;
+  m = 0;
+  last = first + n;
+  for (i = first; i < last; i++) {
+    if (comm_stage == 0) {
+      fp_store[i][0] = buf[m++];
+      fp_store[i][1] = buf[m++];
+      fp_store[i][2] = buf[m++];
+    } else {
+      chi[i] = buf[m++];
+      rho[i] = buf[m++];
+    }
+  }
+}
+
+
+/* ---------------------------------------------------------------------- */
+
+int PairRHEOTension::pack_reverse_comm(int n, int first, double *buf)
+{
+  int i, k, m, last;
+  double **fp_store = compute_interface->fp_store;
+
+  m = 0;
+  last = first + n;
+  for (i = first; i < last; i++) {
+    buf[m++] = fp_store[i][0];
+    buf[m++] = fp_store[i][1];
+    buf[m++] = fp_store[i][2];
+  }
+
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairRHEOTension::unpack_reverse_comm(int n, int *list, double *buf)
+{
+  int i, j, k, m;
+  double **fp_store = compute_interface->fp_store;
+
+  m = 0;
+  for (i = 0; i < n; i++) {
+    j = list[i];
+    fp_store[j][0] += buf[m++];
+    fp_store[j][1] += buf[m++];
+    fp_store[j][2] += buf[m++];
+  }
+}

src/RHEO/pair_rheo_tension.h

@@ -0,0 +1,56 @@
+/* -*- 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 PAIR_CLASS
+// clang-format off
+PairStyle(rheo/tension,PairRHEOTension)
+// clang-format on
+#else
+
+#ifndef LMP_PAIR_RHEO_TENSION_H
+#define LMP_PAIR_RHEO_TENSION_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairRHEOTension : public Pair {
+ public:
+  PairRHEOTension(class LAMMPS *);
+  ~PairRHEOTension() override;
+  void compute(int, int) override;
+  void coeff(int, char **) override;
+  void setup() override;
+  void init_style() override;
+  double init_one(int, int) 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;
+
+ protected:
+  int nmax_store;
+  double **nt, *ct;
+  double *alpha;
+  double hsq, hinv, hinv3;
+
+  void allocate();
+
+  class ComputeRHEOKernel *compute_kernel;
+  class FixRHEO *fix_rheo;
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif