From dfc47a55010a896c99c661f3e13b525f52242eef Mon Sep 17 00:00:00 2001 From: jtclemm Date: Fri, 12 May 2023 23:33:02 -0600 Subject: [PATCH] Fixing various errors, reorganizing status variable --- .../rheo/poiseuille/in.rheo.poiseuille | 2 +- src/.gitignore | 2 + src/RHEO/compute_rheo_grad.cpp | 11 ++-- src/RHEO/compute_rheo_interface.cpp | 10 ++-- src/RHEO/compute_rheo_kernel.cpp | 4 +- src/RHEO/compute_rheo_property_atom.cpp | 8 +-- src/RHEO/compute_rheo_rho_sum.cpp | 1 - src/RHEO/compute_rheo_surface.cpp | 4 +- src/RHEO/compute_rheo_vshift.cpp | 15 ++++-- src/RHEO/fix_rheo.cpp | 34 ++++++------- src/RHEO/fix_rheo.h | 29 ++++++----- src/RHEO/fix_rheo_pressure.cpp | 37 ++------------ src/RHEO/fix_rheo_pressure.h | 3 -- src/RHEO/fix_rheo_thermal.cpp | 22 +++++--- src/RHEO/pair_rheo.cpp | 50 ++++++++++++------- src/RHEO/pair_rheo.h | 1 + 16 files changed, 113 insertions(+), 120 deletions(-) diff --git a/examples/PACKAGES/rheo/poiseuille/in.rheo.poiseuille b/examples/PACKAGES/rheo/poiseuille/in.rheo.poiseuille index af5728c1a3..d0f966c2ce 100644 --- a/examples/PACKAGES/rheo/poiseuille/in.rheo.poiseuille +++ b/examples/PACKAGES/rheo/poiseuille/in.rheo.poiseuille @@ -59,7 +59,7 @@ pair_coeff * * # ------ Fixes & computes ------ # -fix 1 all rheo ${cut} Quintic 0 shift +fix 1 all rheo ${cut} quintic 0 shift fix 2 all rheo/viscosity constant ${eta} fix 3 all rheo/pressure linear fix 4 rig setforce 0.0 0.0 0.0 diff --git a/src/.gitignore b/src/.gitignore index d0fcaf495c..f9794ddb82 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -205,6 +205,8 @@ /compute_rheo_interface.h /compute_rheo_kernel.cpp /compute_rheo_kernel.h +/compute_rheo_property_atom.cpp +/compute_rheo_property_atom.h /compute_rheo_rho_sum.cpp /compute_rheo_rho_sum.h /compute_rheo_surface.cpp diff --git a/src/RHEO/compute_rheo_grad.cpp b/src/RHEO/compute_rheo_grad.cpp index b71fb08d78..92ac108377 100644 --- a/src/RHEO/compute_rheo_grad.cpp +++ b/src/RHEO/compute_rheo_grad.cpp @@ -125,7 +125,7 @@ void ComputeRHEOGrad::init_list(int /*id*/, NeighList *ptr) void ComputeRHEOGrad::compute_peratom() { - int i, j, k, ii, jj, jnum, itype, jtype, a, b; + int i, j, k, ii, jj, jnum, itype, jtype, a, b, fluidi, fluidj; double xtmp, ytmp, ztmp, delx, dely, delz; double rsq, imass, jmass; double rhoi, rhoj, Voli, Volj, drho, dT, deta; @@ -183,6 +183,7 @@ void ComputeRHEOGrad::compute_peratom() vi[1] = v[i][1]; vi[2] = v[i][2]; itype = type[i]; + fluidi = !(status[i] & PHASECHECK); jlist = firstneigh[i]; jnum = numneigh[i]; @@ -197,6 +198,8 @@ void ComputeRHEOGrad::compute_peratom() rsq = delx * delx + dely * dely + delz * delz; if (rsq < cutsq) { + fluidj = !(status[j] & PHASECHECK); + rhoi = rho[i]; rhoj = rho[j]; @@ -206,13 +209,13 @@ void ComputeRHEOGrad::compute_peratom() // Add corrections for walls if (interface_flag) { - if ((status[i] & STATUS_FLUID) && !(status[j] & STATUS_FLUID)) { + if (fluidi && (!fluidj)) { compute_interface->correct_v(vi, vj, i, j); rhoj = compute_interface->correct_rho(j, i); - } else if (!(status[i] & STATUS_FLUID) && (status[j] & STATUS_FLUID)) { + } else if ((!fluidi) && fluidj) { compute_interface->correct_v(vj, vi, j, i); rhoi = compute_interface->correct_rho(i, j); - } else if (!(status[i] & STATUS_FLUID) && !(status[j] & STATUS_FLUID)) { + } else if ((!fluidi) && (!fluidj)) { rhoi = rho0; rhoj = rho0; } diff --git a/src/RHEO/compute_rheo_interface.cpp b/src/RHEO/compute_rheo_interface.cpp index a3624f9663..ea4916087f 100644 --- a/src/RHEO/compute_rheo_interface.cpp +++ b/src/RHEO/compute_rheo_interface.cpp @@ -131,7 +131,7 @@ void ComputeRHEOInterface::compute_peratom() } for (i = 0; i < nall; i++) { - if (!(status[i] & STATUS_FLUID)) rho[i] = 0.0; + if (status[i] & PHASECHECK) rho[i] = 0.0; normwf[i] = 0.0; norm[i] = 0.0; chi[i] = 0.0; @@ -143,7 +143,7 @@ void ComputeRHEOInterface::compute_peratom() ytmp = x[i][1]; ztmp = x[i][2]; itype = type[i]; - fluidi = status[i] & STATUS_FLUID; + fluidi = !(status[i] & PHASECHECK); jlist = firstneigh[i]; jnum = numneigh[i]; @@ -158,7 +158,7 @@ void ComputeRHEOInterface::compute_peratom() if (rsq < cutsq) { jtype = type[j]; - fluidj = status[j] & STATUS_FLUID; + fluidj = !(status[j] & PHASECHECK); w = compute_kernel->calc_w_quintic(i, j, delx, dely, delz, sqrt(rsq)); status_match = 0; @@ -202,7 +202,7 @@ void ComputeRHEOInterface::compute_peratom() if (norm[i] != 0.0) chi[i] /= norm[i]; // Recalculate rho for non-fluid particles - if (!(status[i] & STATUS_FLUID)) { + if (status[i] & PHASECHECK) { if (normwf[i] != 0.0) { // Stores rho for solid particles 1+Pw in Adami Adams 2012 rho[i] = MAX(EPSILON, rho0 + (rho[i] / normwf[i]) * csq_inv); @@ -289,7 +289,7 @@ void ComputeRHEOInterface::unpack_reverse_comm(int n, int *list, double *buf) j = list[i]; norm[j] += buf[m++]; chi[j] += buf[m++]; - if (!(status[j] & STATUS_FLUID)){ + if (status[j] & PHASECHECK){ normwf[j] += buf[m++]; rho[j] += buf[m++]; } else { diff --git a/src/RHEO/compute_rheo_kernel.cpp b/src/RHEO/compute_rheo_kernel.cpp index 09d807d50d..52380a4337 100644 --- a/src/RHEO/compute_rheo_kernel.cpp +++ b/src/RHEO/compute_rheo_kernel.cpp @@ -533,7 +533,7 @@ void ComputeRHEOKernel::compute_peratom() w = calc_w_quintic(i,j,dx[0],dx[1],dx[2],r); rhoj = rho[j]; if (interface_flag) - if (!(status[j] & STATUS_FLUID)) + if (status[j] & PHASECHECK) rhoj = compute_interface->correct_rho(j,i); vj = mass[type[j]] / rhoj; @@ -582,7 +582,7 @@ void ComputeRHEOKernel::compute_peratom() rhoj = rho[j]; if (interface_flag) - if (!(status[j] & STATUS_FLUID)) + if (status[j] & PHASECHECK) rhoj = compute_interface->correct_rho(j,i); vj = mass[type[j]] / rhoj; diff --git a/src/RHEO/compute_rheo_property_atom.cpp b/src/RHEO/compute_rheo_property_atom.cpp index 3cd5d468b2..880ae5d64f 100644 --- a/src/RHEO/compute_rheo_property_atom.cpp +++ b/src/RHEO/compute_rheo_property_atom.cpp @@ -212,10 +212,8 @@ void ComputeRHEOPropertyAtom::pack_phase(int n) int *mask = atom->mask; int nlocal = atom->nlocal; - int inverse_mask = ~PHASEMASK; - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) buf[n] = (status[i] & inverse_mask); + if (mask[i] & groupbit) buf[n] = (status[i] & PHASECHECK); else buf[n] = 0.0; n += nvalues; } @@ -244,10 +242,8 @@ void ComputeRHEOPropertyAtom::pack_surface(int n) int *mask = atom->mask; int nlocal = atom->nlocal; - int inverse_mask = ~SURFACEMASK; - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) buf[n] = (status[i] & inverse_mask); + if (mask[i] & groupbit) buf[n] = (status[i] & SURFACECHECK); else buf[n] = 0.0; n += nvalues; } diff --git a/src/RHEO/compute_rheo_rho_sum.cpp b/src/RHEO/compute_rheo_rho_sum.cpp index 726d876ea1..0a2096a2b9 100644 --- a/src/RHEO/compute_rheo_rho_sum.cpp +++ b/src/RHEO/compute_rheo_rho_sum.cpp @@ -74,7 +74,6 @@ void ComputeRHEORhoSum::compute_peratom() double **x = atom->x; double *rho = atom->rho; int *type = atom->type; - int *status = atom->status; double *mass = atom->mass; int newton = force->newton; diff --git a/src/RHEO/compute_rheo_surface.cpp b/src/RHEO/compute_rheo_surface.cpp index 180c430dd1..f6b93ee551 100644 --- a/src/RHEO/compute_rheo_surface.cpp +++ b/src/RHEO/compute_rheo_surface.cpp @@ -183,7 +183,7 @@ void ComputeRHEOSurface::compute_peratom() jlist = firstneigh[i]; jnum = numneigh[i]; itype = type[i]; - fluidi = status[i] & STATUS_FLUID; + fluidi = !(status[i] & PHASECHECK); for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; @@ -196,7 +196,7 @@ void ComputeRHEOSurface::compute_peratom() rsq = lensq3(dx); if (rsq < cutsq) { jtype = type[j]; - fluidj = status[j] & STATUS_FLUID; + fluidj = !(status[j] & PHASECHECK); rhoi = rho[i]; rhoj = rho[j]; diff --git a/src/RHEO/compute_rheo_vshift.cpp b/src/RHEO/compute_rheo_vshift.cpp index 3d3914436e..0521ff16c0 100644 --- a/src/RHEO/compute_rheo_vshift.cpp +++ b/src/RHEO/compute_rheo_vshift.cpp @@ -91,7 +91,7 @@ void ComputeRHEOVShift::compute_peratom() 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}; + double dx[3], vi[3], vj[3]; int dim = domain->dimension; int *jlist; @@ -123,6 +123,11 @@ void ComputeRHEOVShift::compute_peratom() for (a = 0; a < dim; a++) vshift[i][a] = 0.0; + for (a = 0; a < 3; a++) { + vi[a] = 0.0; + vj[a] = 0.0; + } + for (ii = 0; ii < inum; ii++) { i = ilist[ii]; xtmp = x[i][0]; @@ -132,15 +137,15 @@ void ComputeRHEOVShift::compute_peratom() jlist = firstneigh[i]; jnum = numneigh[i]; imass = mass[itype]; - fluidi = status[i] & STATUS_FLUID; + fluidi = !(status[i] & PHASECHECK); for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; j &= NEIGHMASK; - fluidj = status[j] & STATUS_FLUID; + fluidj = !(status[j] & PHASECHECK); if ((!fluidi) && (!fluidj)) continue; - if (!(status[i] & STATUS_SHIFT) && !(status[j] & STATUS_SHIFT)) continue; + if ((status[i] & STATUS_NO_SHIFT) && (status[j] & STATUS_NO_SHIFT)) continue; dx[0] = xtmp - x[j][0]; dx[1] = ytmp - x[j][1]; @@ -154,7 +159,7 @@ void ComputeRHEOVShift::compute_peratom() r = sqrt(rsq); rinv = 1 / r; - for (a = 0; a < dim; a ++) { + for (a = 0; a < dim; a++) { vi[a] = v[i][a]; vj[a] = v[j][a]; } diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index 04e6c08917..fd436cab6e 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -210,6 +210,8 @@ void FixRHEO::setup_pre_force(int /*vflag*/) // Calculate surfaces if (surface_flag) compute_surface->compute_peratom(); + + pre_force(0); } /* ---------------------------------------------------------------------- */ @@ -258,8 +260,6 @@ void FixRHEO::setup(int /*vflag*/) 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"); - - pre_force(0); } /* ---------------------------------------------------------------------- */ @@ -283,7 +283,8 @@ void FixRHEO::initial_integrate(int /*vflag*/) double **gradr = compute_grad->gradr; double **gradv = compute_grad->gradv; double **vshift; - if (shift_flag) compute_vshift->vshift; + if (shift_flag) + vshift = compute_vshift->vshift; int nlocal = atom->nlocal; int rmass_flag = atom->rmass_flag; @@ -294,7 +295,7 @@ void FixRHEO::initial_integrate(int /*vflag*/) //Density Half-step for (i = 0; i < nlocal; i++) { - if (status[i] & STATUS_NO_FORCE) continue; + if (status[i] & STATUS_NO_INTEGRATION) continue; if (mask[i] & groupbit) { if (rmass_flag) { @@ -331,8 +332,8 @@ void FixRHEO::initial_integrate(int /*vflag*/) if (!rhosum_flag) { for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - if (status[i] & STATUS_NO_FORCE) continue; - if (!(status[i] & STATUS_FLUID)) continue; + if (status[i] & STATUS_NO_INTEGRATION) continue; + if (status[i] & PHASECHECK) continue; divu = 0; for (a = 0; a < dim; a++) { @@ -348,7 +349,7 @@ void FixRHEO::initial_integrate(int /*vflag*/) compute_vshift->correct_surfaces(); // Could this be moved to preforce after the surface fix runs? for (i = 0; i < nlocal; i++) { - if (!(status[i] & STATUS_SHIFT)) continue; + if (status[i] & STATUS_NO_SHIFT) continue; if (mask[i] & groupbit) { for (a = 0; a < dim; a++) { @@ -387,18 +388,13 @@ void FixRHEO::pre_force(int /*vflag*/) if (shift_flag) compute_vshift->compute_peratom(); - // Remove extra shifting/no force options + // Remove temporary options int *mask = atom->mask; int *status = atom->status; int nall = atom->nlocal + atom->nghost; - for (int i = 0; i < nall; i++) { - if (mask[i] & groupbit) { - status[i] &= ~STATUS_NO_FORCE; - - if (status[i] & STATUS_FLUID) - status[i] &= ~STATUS_SHIFT; - } - } + for (int i = 0; i < nall; i++) + if (mask[i] & groupbit) + status[i] &= OPTIONSMASK; // Calculate surfaces, update status if (surface_flag) compute_surface->compute_peratom(); @@ -433,7 +429,7 @@ void FixRHEO::final_integrate() // Update velocity for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - if (status[i] & STATUS_NO_FORCE) continue; + if (status[i] & STATUS_NO_INTEGRATION) continue; if (rmass_flag) { dtfm = dtf / rmass[i]; @@ -451,8 +447,8 @@ void FixRHEO::final_integrate() if (!rhosum_flag) { for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - if (status[i] & STATUS_NO_FORCE) continue; - if (!(status[i] & STATUS_FLUID)) continue; + if (status[i] & STATUS_NO_INTEGRATION) continue; + if (status[i] & PHASECHECK) continue; divu = 0; for (a = 0; a < dim; a++) { diff --git a/src/RHEO/fix_rheo.h b/src/RHEO/fix_rheo.h index 743e418f9a..0dbc8db78b 100644 --- a/src/RHEO/fix_rheo.h +++ b/src/RHEO/fix_rheo.h @@ -74,22 +74,27 @@ namespace RHEO_NS { // Status variables enum Status{ // Phase status - STATUS_FLUID = 1 << 0, - STATUS_SOLID = 1 << 1, - STATUS_REACTIVE = 1 << 2, - STATUS_FREEZING = 1 << 3, + STATUS_SOLID = 1 << 0, + STATUS_REACTIVE = 1 << 1, + // Surface status - STATUS_BULK = 1 << 4, - STATUS_LAYER = 1 << 5, - STATUS_SURFACE = 1 << 6, - STATUS_SPLASH = 1 << 7, + STATUS_BULK = 1 << 2, + STATUS_LAYER = 1 << 3, + STATUS_SURFACE = 1 << 4, + STATUS_SPLASH = 1 << 5, + // Temporary status options - reset in preforce - STATUS_SHIFT = 1 << 8, - STATUS_NO_FORCE = 1 << 9 + STATUS_NO_SHIFT = 1 << 6, + STATUS_NO_INTEGRATION = 1 << 7, + STATUS_FREEZING = 1 << 8 }; - #define PHASEMASK 0xFFFFFFF0; - #define SURFACEMASK 0xFFFFFF0F; + // Masks and their inverses + #define PHASEMASK 0xFFFFFFFC + #define PHASECHECK 0x00000003 + #define SURFACEMASK 0xFFFFFFC3 + #define SURFACECHECK 0x0000003C + #define OPTIONSMASK 0xFFFFFE3F } // namespace RHEO_NS } // namespace LAMMPS_NS diff --git a/src/RHEO/fix_rheo_pressure.cpp b/src/RHEO/fix_rheo_pressure.cpp index 63a6995646..ff206937f4 100644 --- a/src/RHEO/fix_rheo_pressure.cpp +++ b/src/RHEO/fix_rheo_pressure.cpp @@ -43,9 +43,7 @@ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) : if (narg < 4) error->all(FLERR,"Illegal fix command"); pressure_style = NONE; - comm_forward = 1; - nmax_store = 0; // Currently can only have one instance of fix rheo/pressure if (igroup != 0) @@ -73,10 +71,6 @@ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) : FixRHEOPressure::~FixRHEOPressure() { - // Remove custom property if it exists - int tmp1, tmp2, index; - index = atom->find_custom("rheo_pressure", tmp1, tmp2); - if (index != -1) atom->remove_custom(index, 1, 0); } /* ---------------------------------------------------------------------- */ @@ -110,20 +104,6 @@ void FixRHEOPressure::init() void FixRHEOPressure::setup_pre_force(int /*vflag*/) { fix_rheo->pressure_fix_defined = 1; - - // Create pressure array if it doesn'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 - - int tmp1, tmp2; - int index = atom->find_custom("rheo_pressure", tmp1, tmp2); - if (index == -1) { - index = atom->add_custom("rheo_pressure", 1, 0); - nmax_store = atom->nmax; - } - pressure = atom->dvector[index]; - pre_force(0); } @@ -138,14 +118,10 @@ void FixRHEOPressure::pre_force(int /*vflag*/) int *mask = atom->mask; double *rho = atom->rho; + double *pressure = atom->pressure; int nlocal = atom->nlocal; - if (nmax_store < atom->nmax) { - memory->grow(pressure, atom->nmax, "atom:rheo_pressure"); - nmax_store = atom->nmax; - } - for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { if (pressure_style == LINEAR) { @@ -170,6 +146,7 @@ int FixRHEOPressure::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/) { int i,j,k,m; + double *pressure = atom->pressure; m = 0; for (i = 0; i < n; i++) { @@ -184,6 +161,7 @@ int FixRHEOPressure::pack_forward_comm(int n, int *list, double *buf, void FixRHEOPressure::unpack_forward_comm(int n, int first, double *buf) { int i, k, m, last; + double *pressure = atom->pressure; m = 0; last = first + n; @@ -210,12 +188,3 @@ double FixRHEOPressure::calc_pressure(double rho) } return rho; } - -/* ---------------------------------------------------------------------- */ - -double FixRHEOPressure::memory_usage() -{ - double bytes = 0.0; - bytes += (size_t) nmax_store * sizeof(double); - return bytes; -} diff --git a/src/RHEO/fix_rheo_pressure.h b/src/RHEO/fix_rheo_pressure.h index e8f7f3cb88..cbcb495244 100644 --- a/src/RHEO/fix_rheo_pressure.h +++ b/src/RHEO/fix_rheo_pressure.h @@ -34,14 +34,11 @@ class FixRHEOPressure : public Fix { void pre_force(int) override; int pack_forward_comm(int, int *, double *, int, int *) override; void unpack_forward_comm(int, int, double *) override; - double memory_usage() override; double calc_pressure(double); private: double c_cubic, csq, rho0, rho0inv; - double *pressure; int pressure_style; - int nmax_store; class FixRHEO *fix_rheo; }; diff --git a/src/RHEO/fix_rheo_thermal.cpp b/src/RHEO/fix_rheo_thermal.cpp index ec39a13311..de88b4f8d0 100644 --- a/src/RHEO/fix_rheo_thermal.cpp +++ b/src/RHEO/fix_rheo_thermal.cpp @@ -205,7 +205,7 @@ void FixRHEOThermal::initial_integrate(int /*vflag*/) nlocal = atom->nfirst; for (i = 0; i < nlocal; i++) { - if (!(status[i] & STATUS_SHIFT)) continue; + if (status[i] & STATUS_NO_SHIFT) continue; if (mask[i] & groupbit) { for (a = 0; a < dim; a++) { @@ -231,7 +231,7 @@ void FixRHEOThermal::post_integrate() //Integrate temperature and check status for (int i = 0; i < atom->nlocal; i++) { if (mask[i] & groupbit) { - if (status[i] & STATUS_NO_FORCE) continue; + if (status[i] & STATUS_NO_INTEGRATION) continue; cvi = calc_cv(i); temperature[i] += dtf * heatflow[i] / cvi; @@ -245,11 +245,17 @@ void FixRHEOThermal::post_integrate() } if (Ti > Tci) { - status[i] &= PHASEMASK; - status[i] |= STATUS_FLUID; - } else if (!(status[i] & STATUS_SOLID)) { - status[i] &= PHASEMASK; - status[i] |= STATUS_FREEZING; + // If solid, melt + if (status[i] & STATUS_SOLID) { + status[i] &= PHASEMASK; + } + } else { + // If fluid, freeze + if (!(status[i] & STATUS_SOLID)) { + status[i] &= PHASEMASK; + status[i] |= STATUS_SOLID; + status[i] |= STATUS_FREEZING; + } } } } @@ -300,7 +306,7 @@ void FixRHEOThermal::final_integrate() //Integrate temperature and check status for (int i = 0; i < atom->nlocal; i++) { if (mask[i] & groupbit) { - if (status[i] & STATUS_NO_FORCE) continue; + if (status[i] & STATUS_NO_INTEGRATION) continue; cvi = calc_cv(i); temperature[i] += dtf * heatflow[i] / cvi; diff --git a/src/RHEO/pair_rheo.cpp b/src/RHEO/pair_rheo.cpp index 0d041b1e30..0930f28f98 100644 --- a/src/RHEO/pair_rheo.cpp +++ b/src/RHEO/pair_rheo.cpp @@ -120,6 +120,16 @@ void PairRHEO::compute(int eflag, int vflag) numneigh = list->numneigh; firstneigh = list->firstneigh; + for (a = 0; a < 3; a++) { + vi[a] = 0.0; + vj[a] = 0.0; + du[a] = 0.0; + fv[a] = 0.0; + dfp[a] = 0.0; + fsolid[a] = 0.0; + ft[0] = 0.0; + } + // loop over neighbors of my atoms for (ii = 0; ii < inum; ii++) { @@ -132,7 +142,7 @@ void PairRHEO::compute(int eflag, int vflag) jnum = numneigh[i]; imass = mass[itype]; etai = viscosity[i]; - fluidi = status[i] & STATUS_FLUID; + fluidi = !(status[i] & PHASECHECK); if (thermal_flag) { kappai = conductivity[i]; Ti = temperature[i]; @@ -154,7 +164,7 @@ void PairRHEO::compute(int eflag, int vflag) jmass = mass[jtype]; etaj = viscosity[j]; - fluidj = status[j] & STATUS_FLUID; + fluidj = !(status[j] & PHASECHECK); if (thermal_flag) { Tj = temperature[j]; kappaj = conductivity[j]; @@ -186,25 +196,27 @@ void PairRHEO::compute(int eflag, int vflag) Pi = pressure[i]; Pj = pressure[j]; fmag = 0; - if (fluidi && (!fluidj)) { - compute_interface->correct_v(vi, vj, i, j); - rhoj = compute_interface->correct_rho(j, i); - Pj = fix_pressure->calc_pressure(rhoj); + if (interface_flag) { + if (fluidi && (!fluidj)) { + compute_interface->correct_v(vi, vj, i, j); + rhoj = compute_interface->correct_rho(j, i); + Pj = fix_pressure->calc_pressure(rhoj); - if ((chi[j] > 0.9) && (r < (h * 0.5))) - fmag = (chi[j] - 0.9) * (h * 0.5 - r) * rho0 * csq * h * rinv; + 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 = fix_pressure->calc_pressure(rhoi); + } else if ((!fluidi) && fluidj) { + compute_interface->correct_v(vj, vi, j, i); + rhoi = compute_interface->correct_rho(i, j); + Pi = fix_pressure->calc_pressure(rhoi); - if (chi[i] > 0.9 && r < (h * 0.5)) - fmag = (chi[i] - 0.9) * (h * 0.5 - r) * rho0 * csq * h * rinv; + 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 = rho0; - rhoj = rho0; + } else if ((!fluidi) && (!fluidj)) { + rhoi = rho0; + rhoj = rho0; + } } // Repel if close to inner solid particle @@ -234,7 +246,7 @@ void PairRHEO::compute(int eflag, int vflag) sub3(vi, vj, du); //Add artificial viscous pressure if required - if (artificial_visc_flag && pair_avisc_flag){ + if (artificial_visc_flag && pair_avisc_flag) { //Interpolate velocities to midpoint and use this difference for artificial viscosity for (a = 0; a < dim; a++) for (b = 0; b < dim; b++) @@ -328,6 +340,7 @@ void PairRHEO::compute(int eflag, int vflag) } } } + if (vflag_fdotr) virial_fdotr_compute(); } @@ -421,6 +434,7 @@ void PairRHEO::setup() compute_grad = fix_rheo->compute_grad; compute_interface = fix_rheo->compute_interface; thermal_flag = fix_rheo->thermal_flag; + interface_flag = fix_rheo->interface_flag; csq = fix_rheo->csq; rho0 = fix_rheo->rho0; diff --git a/src/RHEO/pair_rheo.h b/src/RHEO/pair_rheo.h index 49aa1ad025..b30b0c3c04 100644 --- a/src/RHEO/pair_rheo.h +++ b/src/RHEO/pair_rheo.h @@ -42,6 +42,7 @@ class PairRHEO : public Pair { int artificial_visc_flag; int rho_damp_flag; int thermal_flag; + int interface_flag; void allocate();