From 009a976ae2bb309bb3ab0da658451ebe201392ab Mon Sep 17 00:00:00 2001 From: jtclemm Date: Sat, 9 Nov 2024 21:47:38 -0700 Subject: [PATCH] Adding multitype correction to rheo vshfit --- src/RHEO/compute_rheo_vshift.cpp | 344 +++++++++++++++++++++++++++++-- src/RHEO/compute_rheo_vshift.h | 10 +- src/RHEO/fix_rheo.cpp | 44 ++-- src/RHEO/fix_rheo.h | 4 +- 4 files changed, 365 insertions(+), 37 deletions(-) diff --git a/src/RHEO/compute_rheo_vshift.cpp b/src/RHEO/compute_rheo_vshift.cpp index 8581d9ea8f..0047521e37 100644 --- a/src/RHEO/compute_rheo_vshift.cpp +++ b/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_interface(nullptr), compute_kernel(nullptr), compute_surface(nullptr) + Compute(lmp, narg, arg), vshift(nullptr), ct(nullptr), wsame(nullptr), cgradt(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; - memory->create(vshift, nmax_store, 3, "rheo:vshift"); + nmax_store = 0; } /* ---------------------------------------------------------------------- */ @@ -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,16 +312,277 @@ 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++) { - buf[m++] = vshift[i][0]; - buf[m++] = vshift[i][1]; - buf[m++] = vshift[i][2]; + 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,14 +591,28 @@ 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; - for (i = 0; i < n; i++) { - j = list[i]; - vshift[j][0] += buf[m++]; - vshift[j][1] += buf[m++]; - vshift[j][2] += buf[m++]; + 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; } diff --git a/src/RHEO/compute_rheo_vshift.h b/src/RHEO/compute_rheo_vshift.h index 485c6525f3..5d763ab008 100644 --- a/src/RHEO/compute_rheo_vshift.h +++ b/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; diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index a978333aba..2241dbc585 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/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(); + compute_kernel->compute_coordination(); // Needed for rho sum - if (rhosum_flag) - compute_rhosum->compute_peratom(); + if (rhosum_flag) compute_rhosum->compute_peratom(); compute_kernel->compute_peratom(); - // Note on first setup, have no forces for pressure to reference - if (interface_flag) + 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(); - // 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) - status[i] &= OPTIONSMASK; + if (mask[i] & groupbit) status[i] &= OPTIONSMASK; // Calculate surfaces, update status if (surface_flag) { compute_surface->compute_peratom(); - if (shift_flag) - compute_vshift->correct_surfaces(); + if (shift_flag) compute_vshift->correct_surfaces(); } } diff --git a/src/RHEO/fix_rheo.h b/src/RHEO/fix_rheo.h index 4b75299f27..bb2d670f33 100644 --- a/src/RHEO/fix_rheo.h +++ b/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;