From 79ddd1445f051fdc7c78bb03123855a94bef2ad9 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Sun, 26 Mar 2023 17:04:56 -0600 Subject: [PATCH] Misc clean ups, initial draft of interface --- src/RHEO/atom_vec_rheo.cpp | 5 + src/RHEO/compute_rheo_grad.cpp | 27 +- src/RHEO/compute_rheo_grad.h | 1 + src/RHEO/compute_rheo_interface.cpp | 432 ++++++++++++++++++++++++++++ src/RHEO/compute_rheo_interface.h | 61 ++++ src/RHEO/compute_rheo_kernel.cpp | 29 +- src/RHEO/compute_rheo_vshift.cpp | 95 +++--- src/RHEO/compute_rheo_vshift.h | 7 +- src/RHEO/fix_rheo.cpp | 9 +- src/RHEO/fix_rheo_pressure.cpp | 11 +- src/RHEO/fix_rheo_thermal.cpp | 11 +- src/RHEO/fix_rheo_viscosity.cpp | 11 +- src/RHEO/pair_rheo.cpp | 5 + 13 files changed, 634 insertions(+), 70 deletions(-) create mode 100644 src/RHEO/compute_rheo_interface.cpp create mode 100644 src/RHEO/compute_rheo_interface.h diff --git a/src/RHEO/atom_vec_rheo.cpp b/src/RHEO/atom_vec_rheo.cpp index 9effad685c..a8496a18e9 100644 --- a/src/RHEO/atom_vec_rheo.cpp +++ b/src/RHEO/atom_vec_rheo.cpp @@ -11,6 +11,11 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: + Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU) +----------------------------------------------------------------------- */ + #include "atom_vec_rheo.h" #include "atom.h" diff --git a/src/RHEO/compute_rheo_grad.cpp b/src/RHEO/compute_rheo_grad.cpp index ce351048be..03e76b194a 100644 --- a/src/RHEO/compute_rheo_grad.cpp +++ b/src/RHEO/compute_rheo_grad.cpp @@ -11,6 +11,11 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: + Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU) +----------------------------------------------------------------------- */ + #include "compute_rheo_grad.h" #include "atom.h" @@ -457,4 +462,24 @@ void ComputeRHEOGrad::grow_arrays(int nmax) if (eta_flag) memory->grow(gradn, nmax, dim, "atom:rheo_grad_eta"); nmax_old = nmax; -} \ No newline at end of file +} + +/* ---------------------------------------------------------------------- */ + +double ComputeRHEOKernel::memory_usage() +{ + double bytes = 0.0; + if (velocity_flag) + bytes = (size_t) nmax_old * dim * dim * sizeof(double); + + if (rho_flag) + bytes = (size_t) nmax_old * dim * sizeof(double); + + if (temperature_flag) + bytes = (size_t) nmax_old * dim * sizeof(double); + + if (eta_flag) + bytes = (size_t) nmax_old * dim * sizeof(double); + + return bytes; +} diff --git a/src/RHEO/compute_rheo_grad.h b/src/RHEO/compute_rheo_grad.h index 220b5813e3..d86efe47ee 100644 --- a/src/RHEO/compute_rheo_grad.h +++ b/src/RHEO/compute_rheo_grad.h @@ -35,6 +35,7 @@ class ComputeRHEOGrad : public Compute { 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 forward_gradients(); void forward_fields(); double **gradv; diff --git a/src/RHEO/compute_rheo_interface.cpp b/src/RHEO/compute_rheo_interface.cpp new file mode 100644 index 0000000000..e76e46d422 --- /dev/null +++ b/src/RHEO/compute_rheo_interface.cpp @@ -0,0 +1,432 @@ +/* ---------------------------------------------------------------------- + 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), Thomas O'Connor (CMU), Eric Palermo (CMU) +----------------------------------------------------------------------- */ + +#include "compute_rheo_interface.h" + +#include "fix_rheo.h" +#include "compute_rheo_kernel.h" +#include "fix_store.h" +#include "fix.h" +#include +#include +#include "atom.h" +#include "update.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; + +#define EPSILON 1e-1 + +/* ---------------------------------------------------------------------- */ + +ComputeRHEOInterface::ComputeRHEOInterface(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), + id_fix_chi(nullptr) +{ + if (narg != 4) error->all(FLERR,"Illegal compute RHEO/chi command"); + + cut = utils::numeric(FLERR,arg[3],false,lmp); + cutsq = cut*cut; + + wall_max = sqrt(3)/12.0*cut; + + nmax = 0; + + comm_forward = 3; + comm_reverse = 4; + + fix_chi = nullptr; + norm = nullptr; + normwf = nullptr; + + // new id = fix-ID + FIX_STORE_ATTRIBUTE + // new fix group = group for this fix + + id_fix_chi = nullptr; + std::string fixcmd = id + std::string("_chi"); + id_fix_chi = new char[fixcmd.size()+1]; + strcpy(id_fix_chi,fixcmd.c_str()); + fixcmd += fmt::format(" all STORE peratom 0 {}", 1); + modify->add_fix(fixcmd); + fix_chi = (FixStore *) modify->fix[modify->nfix-1]; + chi = fix_chi->vstore; +} + +/* ---------------------------------------------------------------------- */ + +ComputeRHEOInterface::~ComputeRHEOInterface() +{ + if (id_fix_chi && modify->nfix) modify->delete_fix(id_fix_chi); + if (modify->nfix) modify->delete_fix("PROPERTY_ATOM_COMP_RHEO_SOLIDS"); + + memory->destroy(norm); + memory->destroy(normwf); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOInterface::init() +{ + // need an occasional full neighbor list + int irequest = neighbor->request(this,instance_me); + neighbor->requests[irequest]->pair = 0; + neighbor->requests[irequest]->compute = 1; + neighbor->requests[irequest]->half = 1; + neighbor->requests[irequest]->full = 0; + //neighbor->requests[irequest]->occasional = 1; //Anticipate needing regulalry + + int icompute = modify->find_compute("rheo_kernel"); + if (icompute == -1) error->all(FLERR, "Using compute/RHEO/chi without compute/RHEO/kernel"); + + compute_kernel = ((ComputeRHEOKernel *) modify->compute[icompute]); + + //Store persistent per atom quantities - need to be exchanged + char **fixarg = new char*[8]; + fixarg[0] = (char *) "PROPERTY_ATOM_COMP_RHEO_SOLIDS"; + fixarg[1] = (char *) "all"; + fixarg[2] = (char *) "property/atom"; + fixarg[3] = (char *) "d_fx"; + fixarg[4] = (char *) "d_fy"; + fixarg[5] = (char *) "d_fz"; + fixarg[6] = (char *) "ghost"; + fixarg[7] = (char *) "yes"; + modify->add_fix(8,fixarg,1); + delete [] fixarg; + + int temp_flag; + index_fx = atom->find_custom("fx", temp_flag); + if ((index_fx < 0) || (temp_flag != 1)) + error->all(FLERR, "Compute rheo/solids can't find fix property/atom fx"); + index_fy = atom->find_custom("fy", temp_flag); + if ((index_fy < 0) || (temp_flag != 1)) + error->all(FLERR, "Compute rheo/solids can't find fix property/atom fy"); + index_fz = atom->find_custom("fz", temp_flag); + if ((index_fz < 0) || (temp_flag != 1)) + error->all(FLERR, "Compute rheo/solids can't find fix property/atom fz"); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOInterface::init_list(int /*id*/, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOInterface::compute_peratom() +{ + int i, j, ii, jj, jnum, itype, jtype, phase_match; + double xtmp, ytmp, ztmp, delx, dely, delz, rsq; + int *jlist; + double w; + + // neighbor list ariables + int inum, *ilist, *numneigh, **firstneigh; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + + double **x = atom->x; + int *type = atom->type; + int newton = force->newton; + int *phase = atom->phase; + double *rho = atom->rho; + double *fx = atom->dvector[index_fx]; + double *fy = atom->dvector[index_fy]; + double *fz = atom->dvector[index_fz]; + double mi, mj; + + //Declare mass pointer to calculate acceleration from force + double *mass = atom->mass; + + cs2 = 1.0; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + if (atom->nmax > nmax) { + nmax = atom->nmax; + fix_chi->grow_arrays(nmax); + chi = fix_chi->vstore; + memory->destroy(norm); + memory->destroy(normwf); + memory->create(norm, nmax, "RHEO/chi:norm"); + memory->create(normwf, nmax, "RHEO/chi:normwf"); + } + + for (i = 0; i < nall; i++) { + if (phase[i] > FixRHEO::FLUID_MAX) rho[i] = 0.0; + normwf[i] = 0.0; + norm[i] = 0.0; + chi[i] = 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]; + mi = mass[type[i]]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx * delx + dely * dely + delz * delz; + + if (rsq < cutsq) { + jtype = type[j]; + mj = mass[type[j]]; + w = compute_kernel->calc_w_quintic(i, j, delx, dely, delz, sqrt(rsq)); + + phase_match = 0; + norm[i] += w; + if ((phase[i] <= FixRHEO::FLUID_MAX && phase[j] <= FixRHEO::FLUID_MAX) + || (phase[i] > FixRHEO::FLUID_MAX && phase[j] > FixRHEO::FLUID_MAX)) { + phase_match = 1; + } + + if (phase_match) { + chi[i] += w; + } else { + if (phase[i] > FixRHEO::FLUID_MAX) { + //speed of sound and rho0 assumed to = 1 (units in density not pressure) + //In general, rho is calculated using the force vector on the wall particle + //fx stores f-fp + rho[i] += w*(cs2*(rho[j] - 1.0) - rho[j]*((-fx[j]/mj+fx[i]/mi)*delx + (-fy[j]/mj+fy[i]/mi)*dely + (-fz[j]/mj+fz[i]/mi)*delz)); + //For the specific taste case whre force on wall particles = 0 + //rho[i] += w*(1.0*(rho[j] - 1.0) + rho[j]*(1e-3*delx)); + normwf[i] += w; + } + } + + if (newton || j < nlocal) { + norm[j] += w; + if (phase_match) { + chi[j] += w; + } else { + if (phase[j] > FixRHEO::FLUID_MAX) { + rho[j] += w*(cs2*(rho[i] - 1.0) + rho[i]*((-fx[i]/mi+fx[j]/mj)*delx + (-fy[i]/mi+fy[j]/mj)*dely + (-fz[i]/mi+fz[j]/mj)*delz)); + normwf[j] += w; + } + } + } + } + } + } + + if (newton) comm->reverse_comm_compute(this); + + for (i = 0; i < nlocal; i++) { + if (norm[i] != 0.0) chi[i] /= norm[i]; + if (normwf[i] != 0.0) { // Only if it's a wall particle + rho[i] = 1.0 + (rho[i] / normwf[i])/cs2; // Stores rho for solid particles 1+Pw in Adami Adams 2012 + if (rho[i] < EPSILON) rho[i] = EPSILON; + } + + if (normwf[i] == 0.0 && phase[i] > FixRHEO::FLUID_MAX) rho[i] = 1.0; + } + + comm_stage = 1; + comm_forward = 2; + comm->forward_comm_compute(this); +} + +/* ---------------------------------------------------------------------- */ + +int ComputeRHEOInterface::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; + double *fx = atom->dvector[index_fx]; + double *fy = atom->dvector[index_fy]; + double *fz = atom->dvector[index_fz]; + + for (i = 0; i < n; i++) { + j = list[i]; + if (comm_stage == 0) { + buf[m++] = fx[j]; + buf[m++] = fy[j]; + buf[m++] = fz[j]; + } else { + buf[m++] = chi[j]; + buf[m++] = rho[j]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOInterface::unpack_forward_comm(int n, int first, double *buf) +{ + int i, k, m, last; + double *rho = atom->rho; + double *fx = atom->dvector[index_fx]; + double *fy = atom->dvector[index_fy]; + double *fz = atom->dvector[index_fz]; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + if (comm_stage == 0) { + fx[i] = buf[m++]; + fy[i] = buf[m++]; + fz[i] = buf[m++]; + } else { + chi[i] = buf[m++]; + rho[i] = buf[m++]; // Won't do anything for fluids + } + } +} + +/* ---------------------------------------------------------------------- */ + +int ComputeRHEOInterface::pack_reverse_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++) { + buf[m++] = norm[i]; + buf[m++] = chi[i]; + buf[m++] = normwf[i]; + buf[m++] = rho[i]; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOInterface::unpack_reverse_comm(int n, int *list, double *buf) +{ + int i,k,j,m; + double *rho = atom->rho; + int *phase = atom->phase; + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + norm[j] += buf[m++]; + chi[j] += buf[m++]; + if (phase[j] > FixRHEO::FLUID_MAX){ + normwf[j] += buf[m++]; + rho[j] += buf[m++]; + } else { + m++; + m++; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOInterface::correct_v(double *vi, double *vj, double *vi_out, int i, int j) +{ + double wall_prefactor, wall_denom, wall_numer; + + wall_numer = 2.0*cut*(chi[i]-0.5); + if (wall_numer < 0) wall_numer = 0; + wall_denom = 2.0*cut*(chi[j]-0.5); + if (wall_denom < wall_max) wall_denom = wall_max; + + wall_prefactor = wall_numer / wall_denom; + + vi_out[0] = (vi[0]-vj[0])*wall_prefactor + vi[0]; + vi_out[1] = (vi[1]-vj[1])*wall_prefactor + vi[1]; + vi_out[2] = (vi[2]-vj[2])*wall_prefactor + vi[2]; +} + +/* ---------------------------------------------------------------------- */ + +double ComputeRHEOInterface::correct_rho(int i, int j) // i is wall, j is fluid +{ + //In future may depend on atom type j's pressure equation + return atom->rho[i]; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOInterface::store_forces() +{ + double *fx = atom->dvector[index_fx]; + double *fy = atom->dvector[index_fy]; + double *fz = atom->dvector[index_fz]; + double **f = atom->f; + double **fp = atom->fp; + int *mask = atom->mask; + + int flag; + int ifix = modify->find_fix_by_style("setforce"); + if (ifix != -1) { + for (int i = 0; i < atom->nlocal; i++) { + if (mask[i] & modify->fix[ifix]->groupbit) { + fx[i] = f[i][0]; + fy[i] = f[i][1]; + fz[i] = f[i][2]; + } else { + fx[i] = f[i][0] - fp[i][0]; + fy[i] = f[i][1] - fp[i][1]; + fz[i] = f[i][2] - fp[i][2]; + } + } + } else { + for (int i = 0; i < atom->nlocal; i++) { + fx[i] = f[i][0] - fp[i][0]; + fy[i] = f[i][1] - fp[i][1]; + fz[i] = f[i][2] - fp[i][2]; + } + } + + //Forward comm forces -note only needed here b/c property atom will forward otherwise + comm_forward = 3; + comm_stage = 0; + comm->forward_comm_compute(this); +} + + +/* ---------------------------------------------------------------------- + memory usage of local atom-based array +------------------------------------------------------------------------- */ + +double ComputeRHEOInterface::memory_usage() +{ + double bytes = nmax * sizeof(double); + return bytes; +} + diff --git a/src/RHEO/compute_rheo_interface.h b/src/RHEO/compute_rheo_interface.h new file mode 100644 index 0000000000..635f190aed --- /dev/null +++ b/src/RHEO/compute_rheo_interface.h @@ -0,0 +1,61 @@ +/* -*- 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/interface,ComputeRHEOInterface) +// clang-format on +#else + +#ifndef LMP_COMPUTE_RHEO_INTERFACE_H +#define LMP_COMPUTE_RHEO_INTERFACE_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeRHEOInterface : public Compute { + public: + ComputeRHEOInterface(class LAMMPS *, int, char **); + ~ComputeRHEOInterface(); + void init(); + void init_list(int, class NeighList *); + void compute_peratom(); + int pack_forward_comm(int, int *, double *, int, int *); + void unpack_forward_comm(int, int, double *); + int pack_reverse_comm(int, int, double *); + void unpack_reverse_comm(int, int *, double *); + void correct_v(double *, double *, double *, int, int); + double correct_rho(int, int); + double memory_usage(); + void store_forces(); + double *chi; + + private: + int nmax; + double cut, cutsq, cs2; + class NeighList *list; + double *norm, *normwf, wall_max; + + class ComputeRHEOKernel *compute_kernel; + char *id_fix_chi; + class FixStore *fix_chi; + + int index_fx, index_fy, index_fz; + int comm_stage; +}; + +} + +#endif +#endif diff --git a/src/RHEO/compute_rheo_kernel.cpp b/src/RHEO/compute_rheo_kernel.cpp index a3f7ed3da5..15b5f9568b 100644 --- a/src/RHEO/compute_rheo_kernel.cpp +++ b/src/RHEO/compute_rheo_kernel.cpp @@ -1,15 +1,20 @@ /* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov + 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. + 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. -------------------------------------------------------------------------- */ + See the README file in the top-level LAMMPS directory. + ------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing authors: + Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU) +----------------------------------------------------------------------- */ #include "compute_rheo_kernel.h" @@ -798,12 +803,12 @@ void ComputeRHEOKernel::unpack_forward_comm(int n, int first, double *buf) double ComputeRHEOKernel::memory_usage() { double bytes = 0.0; - bytes = (size_t) atom->nmax * sizeof(int); + bytes = (size_t) nmax_old * sizeof(int); if (kernel_type == FixRHEO::CRK0) { - bytes += (size_t) atom->nmax * sizeof(double); + bytes += (size_t) nmax_old * sizeof(double); } else if (correction_order > 0) { - bytes += (size_t) atom->nmax * ncor * Mdim * sizeof(double); + bytes += (size_t) nmax_old * ncor * Mdim * sizeof(double); } return bytes; } diff --git a/src/RHEO/compute_rheo_vshift.cpp b/src/RHEO/compute_rheo_vshift.cpp index 42fa1e8c82..2b7c9394d3 100644 --- a/src/RHEO/compute_rheo_vshift.cpp +++ b/src/RHEO/compute_rheo_vshift.cpp @@ -1,35 +1,36 @@ /* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov + 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. + 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. -------------------------------------------------------------------------- */ + See the README file in the top-level LAMMPS directory. + ------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing authors: + Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU) +----------------------------------------------------------------------- */ #include "compute_rheo_vshift.h" -#include "fix_rheo.h" -#include "compute_rheo_solids.h" + +#include "atom.h" +#include "comm.h" +#include "compute_rheo_interface.h" #include "compute_rheo_grad.h" #include "compute_rheo_kernel.h" -#include "fix_rheo_surface.h" -#include -#include -#include "atom.h" -#include "modify.h" #include "domain.h" +#include "error.h" +#include "fix_rheo.h" +#include "force.h" +#include "memory.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; @@ -44,18 +45,28 @@ ComputeRHEOVShift::ComputeRHEOVShift(LAMMPS *lmp, int narg, char **arg) : 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; + // Create vshift 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_old exceeded + + int tmp1, tmp2; + index_vshift = atom->find_custom("rheo_vshift", tmp1, tmp2); + if (index_vshift == -1) { + index_vshift = atom->add_custom("rheo_vshift", 1, 3); + nmax_old = atom->nmax; + } } /* ---------------------------------------------------------------------- */ ComputeRHEOVShift::~ComputeRHEOVShift() { - memory->destroy(vshift); + // Remove custom property if it exists + int tmp1, tmp2, index; + index = atom->find_custom("rheo_vshift", tmp1, tmp2); + if (index != -1) atom->remove_custom(index_vshift, 1, 3); + } /* ---------------------------------------------------------------------- */ @@ -65,10 +76,8 @@ void ComputeRHEOVShift::init() neighbor->add_request(this, NeighConst::REQ_DEFAULT); surface_flag = 0; - if (fix_rheo->surface_flag) { + 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; @@ -110,6 +119,7 @@ void ComputeRHEOVShift::compute_peratom() int *surface = atom->surface; double *rho = atom->rho; double *mass = atom->mass; + double **vshift = atom->darray[index_vshift]; int newton_pair = force->newton_pair; inum = list->inum; @@ -117,12 +127,10 @@ void ComputeRHEOVShift::compute_peratom() 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; - } + + if (nmax_old < atom->nmax) + memory->grow(vshift, atom->nmax, 3, "atom:rheo_vshift"); + nmax_old = atom->nmax; for (i = 0; i < nall; i++) for (a = 0; a < dim; a++) @@ -224,14 +232,17 @@ void ComputeRHEOVShift::correct_surfaces() int nlocal = atom->nlocal; int i, a, b; int dim = domain->dimension; - int *surface = atom->surface; + double **vshift = atom->darray[index_vshift]; + + int tmp1, tmp2; + int index_nsurf = atom->find_custom("rheo_nsurf", tmp1, tmp2); + if (index_nsurf == -1) error->all(FLERR, "Cannot find rheo nsurf"); + double **nsurf = atom->darray[index_nsurf]; - 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) { + if ((surface[i] & FixRHEO::STATUS_SURFACE) || (surface[i] & FixRHEO::STATUS_LAYER)) { nx = nsurf[i][0]; ny = nsurf[i][1]; vx = vshift[i][0]; @@ -257,6 +268,7 @@ void ComputeRHEOVShift::correct_surfaces() int ComputeRHEOVShift::pack_reverse_comm(int n, int first, double *buf) { int i,m,last; + double **vshift = atom->darray[index_vshift]; m = 0; last = first + n; @@ -273,6 +285,7 @@ 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; + double **vshift = atom->darray[index_vshift]; m = 0; for (i = 0; i < n; i++) { @@ -289,6 +302,6 @@ void ComputeRHEOVShift::unpack_reverse_comm(int n, int *list, double *buf) double ComputeRHEOVShift::memory_usage() { - double bytes = 3 * nmax * sizeof(double); + double bytes = 3 * nmax_old * sizeof(double); return bytes; } diff --git a/src/RHEO/compute_rheo_vshift.h b/src/RHEO/compute_rheo_vshift.h index fceb287510..f93cced31c 100644 --- a/src/RHEO/compute_rheo_vshift.h +++ b/src/RHEO/compute_rheo_vshift.h @@ -37,15 +37,12 @@ class ComputeRHEOVShift : public Compute { void correct_surfaces(); private: - int nmax; + int nmax_old; double dtv, cut, cutsq, cutthird; - int surface_flag; - - double **vshift; + int surface_flag, index_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; diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index 6d16bf2782..25a97bd52a 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -1,7 +1,7 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov + 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 @@ -11,6 +11,11 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: + Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU) +----------------------------------------------------------------------- */ + #include "fix_rheo.h" #include "atom.h" diff --git a/src/RHEO/fix_rheo_pressure.cpp b/src/RHEO/fix_rheo_pressure.cpp index 0a92e10767..4fa4ece0f9 100644 --- a/src/RHEO/fix_rheo_pressure.cpp +++ b/src/RHEO/fix_rheo_pressure.cpp @@ -1,7 +1,7 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov + 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 @@ -11,6 +11,11 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: + Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU) +----------------------------------------------------------------------- */ + #include "fix_rheo_pressure.h" #include "atom.h" @@ -210,6 +215,6 @@ double FixRHEOPressure::calculate_p(double rho) double FixRHEOPressure::memory_usage() { double bytes = 0.0; - bytes += (size_t) atom->nmax * sizeof(double); + bytes += (size_t) nmax_old * sizeof(double); return bytes; } diff --git a/src/RHEO/fix_rheo_thermal.cpp b/src/RHEO/fix_rheo_thermal.cpp index c8171f9beb..3ef9fab853 100644 --- a/src/RHEO/fix_rheo_thermal.cpp +++ b/src/RHEO/fix_rheo_thermal.cpp @@ -1,7 +1,7 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov + 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 @@ -11,6 +11,11 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: + Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU) +----------------------------------------------------------------------- */ + #include "fix_rheo_thermal.h" #include "atom.h" @@ -438,6 +443,6 @@ void FixRHEOThermal::unpack_reverse_comm(int n, int *list, double *buf) double FixRHEOThermal::memory_usage() { double bytes = 0.0; - bytes += (size_t) atom->nmax * sizeof(double); + bytes += (size_t) nmax_old * sizeof(double); return bytes; } diff --git a/src/RHEO/fix_rheo_viscosity.cpp b/src/RHEO/fix_rheo_viscosity.cpp index d81a1e9b4f..3f779aa3c6 100644 --- a/src/RHEO/fix_rheo_viscosity.cpp +++ b/src/RHEO/fix_rheo_viscosity.cpp @@ -1,7 +1,7 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov + 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 @@ -11,6 +11,11 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: + Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU) +----------------------------------------------------------------------- */ + #include "fix_rheo_viscosity.h" #include "atom.h" @@ -256,6 +261,6 @@ void FixRHEOViscosity::unpack_forward_comm(int n, int first, double *buf) double FixRHEOViscosity::memory_usage() { double bytes = 0.0; - bytes += (size_t) atom->nmax * sizeof(double); + bytes += (size_t) nmax_old * sizeof(double); return bytes; } diff --git a/src/RHEO/pair_rheo.cpp b/src/RHEO/pair_rheo.cpp index 0fed0297b4..7d5d0f425f 100644 --- a/src/RHEO/pair_rheo.cpp +++ b/src/RHEO/pair_rheo.cpp @@ -11,6 +11,11 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: + Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU) +----------------------------------------------------------------------- */ + #include "pair_rheo.h" #include "atom.h"