Adding multitype correction to rheo vshfit

src/RHEO/compute_rheo_vshift.cpp

@@ -27,6 +27,7 @@
 #include "error.h"
 #include "fix_rheo.h"
 #include "force.h"
+#include "math_extra.h"
 #include "memory.h"
 #include "neigh_list.h"
 #include "neigh_request.h"
@@ -36,20 +37,22 @@
 
 using namespace LAMMPS_NS;
 using namespace RHEO_NS;
+using namespace MathExtra;
 
 /* ---------------------------------------------------------------------- */
 
 ComputeRHEOVShift::ComputeRHEOVShift(LAMMPS *lmp, int narg, char **arg) :
-    Compute(lmp, narg, arg), vshift(nullptr), fix_rheo(nullptr), rho0(nullptr), list(nullptr),
+    Compute(lmp, narg, arg), vshift(nullptr), ct(nullptr), wsame(nullptr), cgradt(nullptr),
-    compute_interface(nullptr), compute_kernel(nullptr), compute_surface(nullptr)
+    fix_rheo(nullptr), rho0(nullptr), list(nullptr), compute_interface(nullptr),
+    compute_kernel(nullptr), compute_surface(nullptr)
 {
   if (narg != 3) error->all(FLERR, "Illegal compute RHEO/VShift command");
 
+  comm_forward = 0;
   comm_reverse = 3;
   surface_flag = 0;
 
-  nmax_store = atom->nmax;
+  nmax_store = 0;
-  memory->create(vshift, nmax_store, 3, "rheo:vshift");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -73,9 +76,19 @@ void ComputeRHEOVShift::init()
   compute_surface = fix_rheo->compute_surface;
 
   rho0 = fix_rheo->rho0;
+  shift_type = fix_rheo->shift_type;
   cut = fix_rheo->cut;
   cutsq = cut * cut;
   cutthird = cut / 3.0;
+
+  multiphase_flag = fix_rheo->shift_multiphase_flag;
+  if (multiphase_flag) {
+    scale = fix_rheo->shift_scale;
+    wmin = fix_rheo->shift_wmin;
+    cmin = fix_rheo->shift_cmin;
+    comm_forward = 1;
+    comm_reverse = 4;
+  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -120,6 +133,13 @@ void ComputeRHEOVShift::compute_peratom()
 
   if (nmax_store < atom->nmax) {
     memory->grow(vshift, atom->nmax, 3, "rheo:vshift");
+
+    if (multiphase_flag) {
+      memory->grow(ct, atom->nmax, "rheo:ct");
+      memory->grow(cgradt, atom->nmax, 3, "rheo:cgradt");
+      memory->grow(wsame, atom->nmax, "rheo:wsame");
+    }
+
     nmax_store = atom->nmax;
   }
 
@@ -224,7 +244,17 @@ void ComputeRHEOVShift::compute_peratom()
     }
   }
 
-  if (newton_pair) comm->reverse_comm(this);
+  comm_stage = 0;
+  if (newton_pair) comm->reverse_comm(this, 3);
+
+  // Zero any excluded types
+
+  for (i = 0; i < nlocal; i++)
+    if (!shift_type[type[i]])
+      for (a = 0; a < dim; a++)
+        vshift[i][a] = 0.0;
+
+  if (multiphase_flag) correct_interfaces();
 }
 
 /* ---------------------------------------------------------------------- */
@@ -246,7 +276,6 @@ void ComputeRHEOVShift::correct_surfaces()
 
       if (status[i] & PHASECHECK) continue;
 
-      //if ((status[i] & STATUS_SURFACE) || (status[i] & STATUS_LAYER)) {
       if (status[i] & STATUS_SURFACE) {
         nx = nsurface[i][0];
         ny = nsurface[i][1];
@@ -283,17 +312,278 @@ void ComputeRHEOVShift::correct_surfaces()
 
 /* ---------------------------------------------------------------------- */
 
-int ComputeRHEOVShift::pack_reverse_comm(int n, int first, double *buf)
+void ComputeRHEOVShift::correct_interfaces()
 {
-  int m, last;
+  int i, j, a, ii, jj, jnum, itype, jtype;
+  int fluidi, fluidj;
+  double xtmp, ytmp, ztmp, rsq, r, rinv;
+  double w, wp, dr, w0, prefactor;
+  double imass, jmass, voli, volj, rhoi, rhoj;
+  double dx[3];
+  int dim = domain->dimension;
+
+  int *jlist;
+  int inum, *ilist, *numneigh, **firstneigh;
+
+  int *type = atom->type;
+  int *status = atom->rheo_status;
+  int *mask = atom->mask;
+  double **x = atom->x;
+  double *rho = atom->rho;
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
+
+  int nlocal = atom->nlocal;
+  int nall = nlocal + atom->nghost;
+  int newton_pair = force->newton_pair;
+
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  size_t nbytes = nmax_store * sizeof(double);
+  memset(&ct[0], 0, nbytes);
+  memset(&wsame[0], 0, nbytes);
+  memset(&cgradt[0][0], 0, 3 * nbytes);
+  double ctmp, *dWij, *dWji;
+
+  // Calculate color gradient
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    itype = type[i];
+    fluidi = !(status[i] & PHASECHECK);
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+    if (rmass)
+      imass = rmass[i];
+    else
+      imass = mass[itype];
+
+    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);
+
+      if (rsq > cutsq) continue;
+
+      fluidj = !(status[j] & PHASECHECK);
+      jtype = type[j];
+      if (rmass)
+        jmass = rmass[j];
+      else
+        jmass = mass[jtype];
+      r = sqrt(rsq);
+
+      rhoi = rho[i];
+      rhoj = rho[j];
+
+      // Add corrections for walls
+      if (interface_flag) {
+        if (fluidi && (!fluidj)) {
+          rhoj = compute_interface->correct_rho(j);
+        } else if ((!fluidi) && fluidj) {
+          rhoi = compute_interface->correct_rho(i);
+        } else if ((!fluidi) && (!fluidj)) {
+          rhoi = rho0[itype];
+          rhoj = rho0[jtype];
+        }
+      }
+
+      voli = imass / rhoi;
+      volj = jmass / rhoj;
+
+      w = compute_kernel->calc_w(i, j, dx[0], dx[1], dx[2], r);
+
+      if (itype != jtype) ctmp = 1;
+      else ctmp = 0;
+
+      ct[i] += ctmp * volj * w;
+      if (newton_pair || j < nlocal)
+        ct[j] += ctmp * voli * w;
+    }
+  }
+
+  comm_stage = 1;
+  if (newton_pair) comm->reverse_comm(this, 1);
+
+  // Calculate color gradient
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    itype = type[i];
+    fluidi = !(status[i] & PHASECHECK);
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+    imass = mass[itype];
+
+    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);
+
+      if (rsq > cutsq) continue;
+
+      fluidj = !(status[j] & PHASECHECK);
+      jtype = type[j];
+      if (rmass)
+        jmass = rmass[j];
+      else
+        jmass = mass[jtype];
+      r = sqrt(rsq);
+
+      rhoi = rho[i];
+      rhoj = rho[j];
+
+      // Add corrections for walls
+      if (interface_flag) {
+        if (fluidi && (!fluidj)) {
+          rhoj = compute_interface->correct_rho(j);
+        } else if ((!fluidi) && fluidj) {
+          rhoi = compute_interface->correct_rho(i);
+        } else if ((!fluidi) && (!fluidj)) {
+          rhoi = rho0[itype];
+          rhoj = rho0[jtype];
+        }
+      }
+
+      voli = imass / rhoi;
+      volj = jmass / rhoj;
+
+      w = compute_kernel->calc_w(i, j, dx[0], dx[1], dx[2], r);
+      dWij = compute_kernel->dWij;
+      dWji = compute_kernel->dWji;
+
+      if (itype != jtype) ctmp = 1;
+      else ctmp = 0;
+
+      for (a = 0; a < dim; a++) {
+        cgradt[i][a] -= ctmp * volj * dWij[a];
+        if (newton_pair || j < nlocal)
+          cgradt[j][a] -= ctmp * voli * dWji[a];
+      }
+
+      if (itype == jtype) {
+        wsame[i] += w * r;
+        if (newton_pair || j < nlocal)
+          wsame[j] += w * r;
+      }
+    }
+  }
+
+  comm_stage = 2;
+  if (newton_pair) comm->reverse_comm(this, 4);
+  comm->forward_comm(this, 1);
+
+  // Correct shifting at fluid-fluid interface
+  // remove normal shifting component for interfacial particles
+  // Based on Yang, Rakhsha, Hu, & Negrut 2022
+
+  double *vs, ntmp[3];
+  double minv, dot;
+
+  for (i = 0; i < nlocal; i++) {
+
+    // If isolated, just don't shift
+    if (wsame[i] < wmin) {
+      for (a = 0; a < dim; a++)
+        vshift[i][a] = 0.0;
+      continue;
+    }
+
+    if (ct[i] < cmin) continue;
+
+    minv = cgradt[i][0] * cgradt[i][0] + cgradt[i][1] * cgradt[i][1];
+    if (dim == 3) minv += cgradt[i][2] * cgradt[i][2];
+
+    if (minv != 0)
+      minv = 1 / sqrt(minv);
+
+    for (a = 0; a < dim; a++)
+      ntmp[a] = cgradt[i][a] * minv;
+
+    vs = vshift[i];
+    dot = ntmp[0] * vs[0] + ntmp[1] * vs[1];
+    if (dim == 3)
+      dot += ntmp[2] * vs[2];
+
+    // To allowing shifting into the bulk
+    // if (dot > 0.0) continue;
+
+    vshift[i][0] -= (1.0 - scale) * ntmp[0] * dot;
+    vshift[i][1] -= (1.0 - scale) * ntmp[1] * dot;
+    if (dim == 3)
+      vshift[i][2] -= (1.0 - scale) * ntmp[2] * dot;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int ComputeRHEOVShift::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/)
+{
+  int i, j, m;
+  m = 0;
+
+  for (i = 0; i < n; i++) {
+    j = list[i];
+    buf[m++] = wsame[j];
+  }
+
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeRHEOVShift::unpack_forward_comm(int n, int first, double *buf)
+{
+  int i, m, last;
 
   m = 0;
   last = first + n;
-  for (int i = first; i < last; i++) {
+  for (i = first; i < last; i++)
+    wsame[i] = buf[m++];
+}
+
+/* ---------------------------------------------------------------------- */
+
+int ComputeRHEOVShift::pack_reverse_comm(int n, int first, double *buf)
+{
+  int i, m, a, last;
+
+  m = 0;
+  last = first + n;
+  if (comm_stage == 0) {
+    for (i = first; i < last; i++) {
       buf[m++] = vshift[i][0];
       buf[m++] = vshift[i][1];
       buf[m++] = vshift[i][2];
     }
+  } else if (comm_stage == 1) {
+    for (i = first; i < last; i++)
+      buf[m++] = ct[i];
+  } else {
+    for (i = first; i < last; i++) {
+      for (a = 0; a < 3; a++)
+        buf[m++] = cgradt[i][a];
+      buf[m++] = wsame[i];
+    }
+  }
   return m;
 }
 
@@ -301,15 +591,29 @@ int ComputeRHEOVShift::pack_reverse_comm(int n, int first, double *buf)
 
 void ComputeRHEOVShift::unpack_reverse_comm(int n, int *list, double *buf)
 {
-  int i, j, m;
+  int i, j, a, m;
 
   m = 0;
+  if (comm_stage == 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++];
     }
+  } else if (comm_stage == 1) {
+    for (i = 0; i < n; i++) {
+      j = list[i];
+      ct[j] += buf[m++];
+    }
+  } else {
+    for (i = 0; i < n; i++) {
+      j = list[i];
+      for (a = 0; a < 3; a++)
+        cgradt[j][a] += buf[m++];
+      wsame[j] += buf[m++];
+    }
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -319,5 +623,9 @@ void ComputeRHEOVShift::unpack_reverse_comm(int n, int *list, double *buf)
 double ComputeRHEOVShift::memory_usage()
 {
   double bytes = 3 * nmax_store * sizeof(double);
+
+  if (multiphase_flag)
+    bytes += 5 * nmax_store * sizeof(double);
+
   return bytes;
 }

src/RHEO/compute_rheo_vshift.h

@@ -31,19 +31,25 @@ class ComputeRHEOVShift : public Compute {
   void init() override;
   void init_list(int, class NeighList *) override;
   void compute_peratom() override;
+  int pack_forward_comm(int, int *, double *, int, int *) override;
+  void unpack_forward_comm(int, int, double *) override;
   int pack_reverse_comm(int, int, double *) override;
   void unpack_reverse_comm(int, int *, double *) override;
   double memory_usage() override;
   void correct_surfaces();
+  void correct_interfaces();
   double **vshift;
 
   class FixRHEO *fix_rheo;
 
  private:
-  int nmax_store;
+  int nmax_store, comm_stage;
   double dtv, cut, cutsq, cutthird;
-  int surface_flag, interface_flag;
+  double scale, wmin, cmin;
+  int surface_flag, interface_flag, multiphase_flag;
   double *rho0;
+  double *wsame, *ct, **cgradt;
+  int *shift_type;
 
   class NeighList *list;
   class ComputeRHEOInterface *compute_interface;

src/RHEO/fix_rheo.cpp

@@ -70,15 +70,17 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
   surface_flag = 0;
   oxidation_flag = 0;
   self_mass_flag = 0;
-  coordination_flag = 0;
+  shift_multiphase_flag = 0;
 
   int i;
   int n = atom->ntypes;
   memory->create(rho0, n + 1, "rheo:rho0");
   memory->create(csq, n + 1, "rheo:csq");
+  memory->create(shift_type, n + 1, "rheo:shift_type");
   for (i = 1; i <= n; i++) {
     rho0[i] = 1.0;
     csq[i] = 1.0;
+    shift_type[i] = 1;
   }
 
   if (igroup != 0) error->all(FLERR, "fix rheo command requires group all");
@@ -112,6 +114,19 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
   while (iarg < narg) {
     if (strcmp(arg[iarg], "shift") == 0) {
       shift_flag = 1;
+    } else if (strcmp(arg[iarg], "shift/multiphase/scale") == 0) {
+      if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift/multiphase/scale", error);
+      shift_multiphase_flag = 1;
+      shift_scale = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
+      shift_wmin = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
+      shift_cmin = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
+      iarg += 3;
+    } else if (strcmp(arg[iarg], "shift/type") == 0) {
+      if (iarg + n >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift/type", error);
+      for (i = 1; i <= n; i++) {
+        shift_type[i] = utils::logical(FLERR, arg[iarg + i], false, lmp);
+      }
+      iarg += n;
     } else if (strcmp(arg[iarg], "thermal") == 0) {
       thermal_flag = 1;
     } else if (strcmp(arg[iarg], "surface/detection") == 0) {
@@ -128,7 +143,6 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
       } else {
         error->all(FLERR, "Illegal surface/detection option in fix rheo, {}", arg[iarg + 1]);
       }
-
       iarg += 3;
     } else if (strcmp(arg[iarg], "interface/reconstruct") == 0) {
       interface_flag = 1;
@@ -153,6 +167,9 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
     iarg += 1;
   }
 
+  if ((!shift_flag) && shift_multiphase_flag)
+    error->all(FLERR, "Cannot use shift/multiphase/scale without shifting");
+
   if (self_mass_flag && (!rhosum_flag))
     error->all(FLERR, "Cannot use self/mass setting without rho/sum");
 
@@ -174,6 +191,7 @@ FixRHEO::~FixRHEO()
 
   memory->destroy(csq);
   memory->destroy(rho0);
+  memory->destroy(shift_type);
 }
 
 /* ----------------------------------------------------------------------
@@ -376,7 +394,6 @@ void FixRHEO::initial_integrate(int /*vflag*/)
   // Shifting atoms
   if (shift_flag) {
     for (i = 0; i < nlocal; i++) {
-
       if (status[i] & STATUS_NO_SHIFT) continue;
       if (status[i] & PHASECHECK) continue;
 
@@ -399,39 +416,34 @@ void FixRHEO::initial_integrate(int /*vflag*/)
 
 void FixRHEO::pre_force(int /*vflag*/)
 {
-  if (coordination_flag)
+  compute_kernel->compute_coordination();    // Needed for rho sum
-    compute_kernel->compute_coordination();
 
-  if (rhosum_flag)
+  if (rhosum_flag) compute_rhosum->compute_peratom();
-    compute_rhosum->compute_peratom();
 
   compute_kernel->compute_peratom();
 
+  if (interface_flag) {
     // Note on first setup, have no forces for pressure to reference
-  if (interface_flag)
     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();
 
-  // Depends on NO_SHIFT status
+  if (shift_flag) compute_vshift->compute_peratom();
-  if (shift_flag)
-    compute_vshift->compute_peratom();
 
   // Remove temporary options
   int *mask = atom->mask;
   int *status = atom->rheo_status;
   int nall = atom->nlocal + atom->nghost;
   for (int i = 0; i < nall; i++)
-    if (mask[i] & groupbit)
+    if (mask[i] & groupbit) status[i] &= OPTIONSMASK;
-      status[i] &= OPTIONSMASK;
 
   // Calculate surfaces, update status
   if (surface_flag) {
     compute_surface->compute_peratom();
-    if (shift_flag)
+    if (shift_flag) compute_vshift->correct_surfaces();
-      compute_vshift->correct_surfaces();
   }
 }
 

src/RHEO/fix_rheo.h

@@ -41,10 +41,12 @@ class FixRHEO : public Fix {
   // Model parameters
   double cut;
   double *rho0, *csq;
+  int *shift_type;
   int self_mass_flag;
   int zmin_kernel, zmin_surface, zmin_splash;
   int kernel_style, surface_style;
   double divr_surface;
+  double shift_scale, shift_wmin, shift_cmin;
 
   // Accessory fixes/computes
   int thermal_flag;
@@ -53,7 +55,7 @@ class FixRHEO : public Fix {
   int interface_flag;
   int surface_flag;
   int oxidation_flag;
-  int coordination_flag;
+  int shift_multiphase_flag;
 
   int viscosity_fix_defined;
   int pressure_fix_defined;