diff --git a/src/RHEO/atom_vec_rheo.cpp b/src/RHEO/atom_vec_rheo.cpp index de2e2f77ad..ea9e2a3c10 100644 --- a/src/RHEO/atom_vec_rheo.cpp +++ b/src/RHEO/atom_vec_rheo.cpp @@ -33,15 +33,17 @@ AtomVecRHEO::AtomVecRHEO(LAMMPS *lmp) : AtomVec(lmp) forceclearflag = 1; atom->status_flag = 1; + atom->pressure_flag = 1; atom->rho_flag = 1; + atom->viscosity_flag = 1; // strings with peratom variables to include in each AtomVec method // strings cannot contain fields in corresponding AtomVec default strings // order of fields in a string does not matter // except: fields_data_atom & fields_data_vel must match data file - fields_grow = {"status", "rho", "drho"}; - fields_copy = {"status", "rho", "drho"}; + fields_grow = {"status", "rho", "drho", "pressure", "viscosity"}; + fields_copy = {"status", "rho", "drho", "pressure", "viscosity"}; fields_comm = {"status", "rho"}; fields_comm_vel = {"status", "rho"}; fields_reverse = {"drho"}; @@ -49,7 +51,7 @@ AtomVecRHEO::AtomVecRHEO(LAMMPS *lmp) : AtomVec(lmp) fields_border_vel = {"status", "rho"}; fields_exchange = {"status", "rho"}; fields_restart = {"status", "rho"}; - fields_create = {"status", "rho", "drho"}; + fields_create = {"status", "rho", "drho", "pressure", "viscosity"}; fields_data_atom = {"id", "type", "status", "rho", "x"}; fields_data_vel = {"id", "v"}; @@ -64,8 +66,10 @@ AtomVecRHEO::AtomVecRHEO(LAMMPS *lmp) : AtomVec(lmp) void AtomVecRHEO::grow_pointers() { status = atom->status; + pressure = atom->pressure; rho = atom->rho; drho = atom->drho; + viscosity = atom->viscosity; } /* ---------------------------------------------------------------------- @@ -86,6 +90,8 @@ void AtomVecRHEO::force_clear(int n, size_t nbytes) void AtomVecRHEO::data_atom_post(int ilocal) { drho[ilocal] = 0.0; + pressure[ilocal] = 0.0; + viscosity[ilocal] = 0.0; } /* ---------------------------------------------------------------------- @@ -96,8 +102,10 @@ void AtomVecRHEO::data_atom_post(int ilocal) int AtomVecRHEO::property_atom(const std::string &name) { if (name == "status") return 0; - if (name == "rho") return 1; - if (name == "drho") return 2; + if (name == "pressure") return 1; + if (name == "rho") return 2; + if (name == "drho") return 3; + if (name == "viscosity") return 4; return -1; } @@ -123,12 +131,20 @@ void AtomVecRHEO::pack_property_atom(int index, double *buf, int nvalues, int gr } else if (index == 1) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) - buf[n] = rho[i]; + buf[n] = pressure[i]; else buf[n] = 0.0; n += nvalues; } } else if (index == 2) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = rho[i]; + else + buf[n] = 0.0; + n += nvalues; + } + } else if (index == 3) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) buf[n] = drho[i]; @@ -136,5 +152,13 @@ void AtomVecRHEO::pack_property_atom(int index, double *buf, int nvalues, int gr buf[n] = 0.0; n += nvalues; } + } else if (index == 4) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = viscosity[i]; + else + buf[n] = 0.0; + n += nvalues; + } } } diff --git a/src/RHEO/atom_vec_rheo.h b/src/RHEO/atom_vec_rheo.h index bdd617a01d..68cc224ba5 100644 --- a/src/RHEO/atom_vec_rheo.h +++ b/src/RHEO/atom_vec_rheo.h @@ -36,7 +36,7 @@ class AtomVecRHEO : virtual public AtomVec { private: int *status; - double *rho, *drho; + double *pressure, *rho, *drho, *viscosity; }; } // namespace LAMMPS_NS diff --git a/src/RHEO/atom_vec_rheo_thermal.cpp b/src/RHEO/atom_vec_rheo_thermal.cpp new file mode 100644 index 0000000000..de0c7fa5d7 --- /dev/null +++ b/src/RHEO/atom_vec_rheo_thermal.cpp @@ -0,0 +1,200 @@ +/* ---------------------------------------------------------------------- + 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 "atom_vec_rheo_thermal.h" + +#include "atom.h" + +#include + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +AtomVecRHEOThermal::AtomVecRHEOThermal(LAMMPS *lmp) : AtomVec(lmp) +{ + molecular = Atom::ATOMIC; + mass_type = PER_TYPE; + forceclearflag = 1; + + atom->status_flag = 1; + atom->conductivity_flag = 1; + atom->temperature_flag = 1; + atom->heatflow_flag = 1; + atom->pressure_flag = 1; + atom->rho_flag = 1; + atom->viscosity_flag = 1; + + // strings with peratom variables to include in each AtomVec method + // strings cannot contain fields in corresponding AtomVec default strings + // order of fields in a string does not matter + // except: fields_data_atom & fields_data_vel must match data file + + fields_grow = {"status", "rho", "drho", "temperature", "heatflow", "conductivity", "pressure", "viscosity"}; + fields_copy = {"status", "rho", "drho", "temperature", "heatflow", "conductivity", "pressure", "viscosity"}; + fields_comm = {"status", "rho", "temperature"}; + fields_comm_vel = {"status", "rho", "temperature"}; + fields_reverse = {"drho", "heatflow"}; + fields_border = {"status", "rho", "temperature"}; + fields_border_vel = {"status", "rho", "temperature"}; + fields_exchange = {"status", "rho", "temperature"}; + fields_restart = {"status", "rho", "temperature"}; + fields_create = {"status", "rho", "drho", "temperature", "heatflow", "conductivity", "pressure", "viscosity"}; + fields_data_atom = {"id", "type", "status", "rho", "temperature", "x"}; + fields_data_vel = {"id", "v"}; + + setup_fields(); +} + +/* ---------------------------------------------------------------------- + set local copies of all grow ptrs used by this class, except defaults + needed in replicate when 2 atom classes exist and it calls pack_restart() +------------------------------------------------------------------------- */ + +void AtomVecRHEOThermal::grow_pointers() +{ + status = atom->status; + conductivity = atom->conductivity; + temperature = atom->temperature; + heatflow = atom->heatflow; + pressure = atom->pressure; + rho = atom->rho; + drho = atom->drho; + viscosity = atom->viscosity; +} + +/* ---------------------------------------------------------------------- + clear extra forces starting at atom N + nbytes = # of bytes to clear for a per-atom vector +------------------------------------------------------------------------- */ + +void AtomVecRHEOThermal::force_clear(int n, size_t nbytes) +{ + memset(&drho[n], 0, nbytes); + memset(&heatflow[n], 0, nbytes); +} + +/* ---------------------------------------------------------------------- + modify what AtomVec::data_atom() just unpacked + or initialize other atom quantities +------------------------------------------------------------------------- */ + +void AtomVecRHEOThermal::data_atom_post(int ilocal) +{ + drho[ilocal] = 0.0; + heatflow[ilocal] = 0.0; + pressure[ilocal] = 0.0; + viscosity[ilocal] = 0.0; + conductivity[ilocal] = 0.0; +} + +/* ---------------------------------------------------------------------- + assign an index to named atom property and return index + return -1 if name is unknown to this atom style +------------------------------------------------------------------------- */ + +int AtomVecRHEOThermal::property_atom(const std::string &name) +{ + if (name == "status") return 0; + if (name == "rho") return 1; + if (name == "drho") return 2; + if (name == "temperature") return 3; + if (name == "heatflow") return 4; + if (name == "conductivity") return 5; + if (name == "pressure") return 6; + if (name == "viscosity") return 7; + return -1; +} + +/* ---------------------------------------------------------------------- + pack per-atom data into buf for ComputePropertyAtom + index maps to data specific to this atom style +------------------------------------------------------------------------- */ + +void AtomVecRHEOThermal::pack_property_atom(int index, double *buf, int nvalues, int groupbit) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + int n = 0; + + if (index == 0) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = status[i]; + else + buf[n] = 0.0; + n += nvalues; + } + } else if (index == 1) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = rho[i]; + else + buf[n] = 0.0; + n += nvalues; + } + } else if (index == 2) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = drho[i]; + else + buf[n] = 0.0; + n += nvalues; + } + } else if (index == 3) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = temperature[i]; + else + buf[n] = 0.0; + n += nvalues; + } + } else if (index == 4) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = heatflow[i]; + else + buf[n] = 0.0; + n += nvalues; + } + } else if (index == 5) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = conductivity[i]; + else + buf[n] = 0.0; + n += nvalues; + } + } else if (index == 6) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = pressure[i]; + else + buf[n] = 0.0; + n += nvalues; + } + } else if (index == 7) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = viscosity[i]; + else + buf[n] = 0.0; + n += nvalues; + } + } +} diff --git a/src/RHEO/atom_vec_rheo_thermal.h b/src/RHEO/atom_vec_rheo_thermal.h new file mode 100644 index 0000000000..27c6c3c9b5 --- /dev/null +++ b/src/RHEO/atom_vec_rheo_thermal.h @@ -0,0 +1,46 @@ +/* -*- 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 ATOM_CLASS +// clang-format off +AtomStyle(rheo/thermal,AtomVecRHEOThermal); +// clang-format on +#else + +#ifndef LMP_ATOM_VEC_RHEO_THERMAL_H +#define LMP_ATOM_VEC_RHEO_THERMAL_H + +#include "atom_vec.h" + +namespace LAMMPS_NS { + +class AtomVecRHEOThermal : virtual public AtomVec { + public: + AtomVecRHEOThermal(class LAMMPS *); + + void grow_pointers() override; + void force_clear(int, size_t) override; + void data_atom_post(int) override; + int property_atom(const std::string &) override; + void pack_property_atom(int, double *, int, int) override; + + private: + int *status; + double *conductivity, *temperature, *heatflow; + double *pressure, *rho, *drho, *viscosity; +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/RHEO/compute_rheo_grad.cpp b/src/RHEO/compute_rheo_grad.cpp index 606cec7dfc..b71fb08d78 100644 --- a/src/RHEO/compute_rheo_grad.cpp +++ b/src/RHEO/compute_rheo_grad.cpp @@ -84,69 +84,34 @@ ComputeRHEOGrad::ComputeRHEOGrad(LAMMPS *lmp, int narg, char **arg) : } comm_forward = ncomm_grad; + + nmax_store = 0; + grow_arrays(atom->nmax); + } /* ---------------------------------------------------------------------- */ ComputeRHEOGrad::~ComputeRHEOGrad() { - int dim = domain->dimension; - int tmp1, tmp2, index; - - index = atom->find_custom("rheo_grad_v", tmp1, tmp2); - if (index != 1) atom->remove_custom(index, 1, dim * dim); - index = atom->find_custom("rheo_grad_rho", tmp1, tmp2); - if (index != 1) atom->remove_custom(index, 1, dim); - index = atom->find_custom("rheo_grad_t", tmp1, tmp2); - if (index != 1) atom->remove_custom(index, 1, dim); - index = atom->find_custom("rheo_grad_eta", tmp1, tmp2); - if (index != 1) atom->remove_custom(index, 1, dim); + memory->destroy(gradv); + memory->destroy(gradr); + memory->destroy(gradt); + memory->destroy(gradn); } /* ---------------------------------------------------------------------- */ void ComputeRHEOGrad::init() { - neighbor->add_request(this, NeighConst::REQ_DEFAULT); - cut = fix_rheo->cut; cutsq = cut * cut; rho0 = fix_rheo->rho0; + interface_flag = fix_rheo->interface_flag; compute_kernel = fix_rheo->compute_kernel; compute_interface = fix_rheo->compute_interface; - int tmp1, tmp2; - index_visc = atom->find_custom("rheo_viscosity", tmp1, tmp2); - - // Create coordination 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 index; - int dim = domain->dimension; - if (velocity_flag) { - index = atom->add_custom("rheo_grad_v", 1, dim * dim); - gradv = atom->darray[index]; - } - - if (rho_flag) { - index = atom->add_custom("rheo_grad_rho", 1, dim); - gradr = atom->darray[index]; - } - - if (temperature_flag) { - index= atom->add_custom("rheo_grad_temp", 1, dim); - gradt = atom->darray[index]; - } - - if (eta_flag) { - index = atom->add_custom("rheo_grad_eta", 1, dim); - gradn = atom->darray[index]; - } - - nmax_store = 0; - grow_arrays(atom->nmax); + neighbor->add_request(this, NeighConst::REQ_DEFAULT); } /* ---------------------------------------------------------------------- */ @@ -175,7 +140,7 @@ void ComputeRHEOGrad::compute_peratom() double **v = atom->v; double *rho = atom->rho; double *temperature = atom->temperature; - double *viscosity = atom->dvector[index_visc]; + double *viscosity = atom->viscosity; int *status = atom->status; int *type = atom->type; double *mass = atom->mass; @@ -240,15 +205,17 @@ void ComputeRHEOGrad::compute_peratom() vj[2] = v[j][2]; // Add corrections for walls - if ((status[i] & STATUS_FLUID) && !(status[j] & STATUS_FLUID)) { - 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)) { - 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)) { - rhoi = rho0; - rhoj = rho0; + if (interface_flag) { + if ((status[i] & STATUS_FLUID) && !(status[j] & STATUS_FLUID)) { + 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)) { + 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)) { + rhoi = rho0; + rhoj = rho0; + } } Voli = mass[itype] / rhoi; @@ -481,16 +448,16 @@ void ComputeRHEOGrad::grow_arrays(int nmax) { int dim = domain->dimension; if (velocity_flag) - memory->grow(gradv, nmax, dim * dim, "atom:rheo_grad_v"); + memory->grow(gradv, nmax, dim * dim, "rheo:grad_v"); if (rho_flag) - memory->grow(gradr, nmax, dim, "atom:rheo_grad_rho"); + memory->grow(gradr, nmax, dim, "rheo:grad_rho"); if (temperature_flag) - memory->grow(gradt, nmax, dim, "atom:rheo_grad_temp"); + memory->grow(gradt, nmax, dim, "rheo:grad_temp"); if (eta_flag) - memory->grow(gradn, nmax, dim, "atom:rheo_grad_eta"); + memory->grow(gradn, nmax, dim, "rheo:grad_eta"); nmax_store = nmax; } diff --git a/src/RHEO/compute_rheo_grad.h b/src/RHEO/compute_rheo_grad.h index 8c7962a978..af4fecdcfb 100644 --- a/src/RHEO/compute_rheo_grad.h +++ b/src/RHEO/compute_rheo_grad.h @@ -46,14 +46,14 @@ class ComputeRHEOGrad : public Compute { private: int comm_stage, ncomm_grad, ncomm_field, nmax_store; - int index_visc; double cut, cutsq, rho0; - class NeighList *list; + + int velocity_flag, temperature_flag, rho_flag, eta_flag; + int interface_flag; class ComputeRHEOKernel *compute_kernel; class ComputeRHEOInterface *compute_interface; - - int velocity_flag, temperature_flag, rho_flag, eta_flag; + class NeighList *list; void grow_arrays(int); }; diff --git a/src/RHEO/compute_rheo_interface.cpp b/src/RHEO/compute_rheo_interface.cpp index 72dcdc17d7..a3624f9663 100644 --- a/src/RHEO/compute_rheo_interface.cpp +++ b/src/RHEO/compute_rheo_interface.cpp @@ -46,24 +46,35 @@ ComputeRHEOInterface::ComputeRHEOInterface(LAMMPS *lmp, int narg, char **arg) : { if (narg != 3) error->all(FLERR,"Illegal compute rheo/interface command"); - nmax_store = 0; - comm_forward = 3; comm_reverse = 4; + + nmax_store = atom->nmax; + memory->create(chi, nmax_store, "rheo:chi"); + memory->create(norm, nmax_store, "rheo/interface:norm"); + memory->create(normwf, nmax_store, "rheo/interface:normwf"); + + // For fp_store, create an instance of fix property atom + // Need restarts + exchanging with neighbors since it needs to persist + // between timesteps (fix property atom will handle callbacks) + + int tmp1, tmp2; + int index = atom->find_custom("fp_store", tmp1, tmp2); + if (index == -1) { + id_fix_pa = utils::strdup(id + std::string("_fix_property_atom")); + modify->add_fix(fmt::format("{} all property/atom d2_fp_store 3", id_fix_pa)); + index = atom->find_custom("fp_store", tmp1, tmp2); + } + fp_store = atom->darray[index]; } /* ---------------------------------------------------------------------- */ ComputeRHEOInterface::~ComputeRHEOInterface() { - // Remove custom property if it exists - int tmp1, tmp2, index; - index = atom->find_custom("rheo_chi", tmp1, tmp2); - if (index != -1) atom->remove_custom(index, 1, 0); - if (id_fix_pa && modify->nfix) modify->delete_fix(id_fix_pa); delete[] id_fix_pa; - + memory->destroy(chi); memory->destroy(norm); memory->destroy(normwf); } @@ -80,37 +91,6 @@ void ComputeRHEOInterface::init() cutsq = cut * cut; wall_max = sqrt(3.0) / 12.0 * cut; - // Create chi 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 nmax = atom->nmax; - int index = atom->find_custom("rheo_chi", tmp1, tmp2); - if (index == -1) { - index = atom->add_custom("rheo_chi", 1, 0); - memory->destroy(norm); - memory->destroy(normwf); - memory->create(norm, nmax, "rheo/interface:norm"); - memory->create(normwf, nmax, "rheo/interface:normwf"); - nmax_store = nmax; - } - chi = atom->dvector[index]; - - // For fp_store, go ahead and create an instance of fix property atom - // Need restarts + exchanging with neighbors since it needs to persist - // between timesteps (fix property atom will handle callbacks) - - index = atom->find_custom("fp_store", tmp1, tmp2); - if (index == -1) { - id_fix_pa = utils::strdup(id + std::string("_fix_property_atom")); - modify->add_fix(fmt::format("{} all property/atom d2_fp_store 3", id_fix_pa)); - index = atom->find_custom("fp_store", tmp1, tmp2); - } - fp_store = atom->darray[index]; - - // need an occasional half neighbor list neighbor->add_request(this, NeighConst::REQ_DEFAULT); } @@ -145,11 +125,9 @@ void ComputeRHEOInterface::compute_peratom() if (atom->nmax > nmax_store) { nmax_store = atom->nmax; - memory->destroy(norm); - memory->destroy(normwf); - memory->create(norm, nmax_store, "rheo/interface:norm"); - memory->create(normwf, nmax_store, "rheo/interface:normwf"); - memory->grow(chi, nmax_store, "rheo/interface:chi"); + memory->grow(norm, nmax_store, "rheo/interface:norm"); + memory->grow(normwf, nmax_store, "rheo/interface:normwf"); + memory->grow(chi, nmax_store, "rheo:chi"); } for (i = 0; i < nall; i++) { diff --git a/src/RHEO/compute_rheo_kernel.cpp b/src/RHEO/compute_rheo_kernel.cpp index b7ad1c9b3f..9cfa86df7e 100644 --- a/src/RHEO/compute_rheo_kernel.cpp +++ b/src/RHEO/compute_rheo_kernel.cpp @@ -68,7 +68,6 @@ ComputeRHEOKernel::ComputeRHEOKernel(LAMMPS *lmp, int narg, char **arg) : correction_order = 2; } - solid_flag = 0; dim = domain->dimension; comm_forward = 1; @@ -93,11 +92,7 @@ ComputeRHEOKernel::ComputeRHEOKernel(LAMMPS *lmp, int narg, char **arg) : ComputeRHEOKernel::~ComputeRHEOKernel() { - // Remove custom property if it exists - int tmp1, tmp2, index; - index = atom->find_custom("rheo_coordination", tmp1, tmp2); - if (index != -1) atom->remove_custom(index, 0, 0); - + memory->destroy(coordination); memory->destroy(C); memory->destroy(C0); } @@ -112,11 +107,8 @@ void ComputeRHEOKernel::init() if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use compute rheo/kernel"); fix_rheo = dynamic_cast(fixes[0]); - int icompute = modify->find_compute("rheo_interface"); - if (icompute != -1) { - compute_interface = ((ComputeRHEOInterface *) modify->compute[icompute]); - solid_flag = 1; - } + interface_flag = fix_rheo->interface_flag; + compute_interface = fix_rheo->compute_interface; zmin = fix_rheo->zmin_kernel; h = fix_rheo->h; @@ -133,22 +125,8 @@ void ComputeRHEOKernel::init() pre_wp = pre_w * 3.0 * hinv; } - // Create coordination 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 nmax = atom->nmax; - int index = atom->find_custom("rheo_coordination", tmp1, tmp2); - if (index == -1) { - index = atom->add_custom("rheo_coordination", 0, 0); - nmax_store = nmax; - } - coordination = atom->ivector[index]; - - // Create local arrays for kernel arrays, I can't foresee a reason to print - + nmax_store = atom->nmax; + memory->create(coordination, nmax_store, "rheo:coordination"); if (kernel_style == CRK0) { memory->create(C0, nmax_store, "rheo/kernel:C0"); } else if (kernel_style == CRK1) { @@ -499,7 +477,7 @@ void ComputeRHEOKernel::compute_peratom() if (kernel_style == QUINTIC) return; int i, j, ii, jj, inum, jnum, itype, g, a, b, gsl_error; - double xtmp, ytmp, ztmp, r, rsq, w, vj; + double xtmp, ytmp, ztmp, r, rsq, w, vj, rhoj; double dx[3]; gsl_matrix_view gM; @@ -549,10 +527,12 @@ void ComputeRHEOKernel::compute_peratom() if (rsq < hsq) { r = sqrt(rsq); w = calc_w_quintic(i,j,dx[0],dx[1],dx[2],r); - if (!(status[j] & STATUS_FLUID) && solid_flag) { - vj = mass[type[j]] / compute_interface->correct_rho(j,i); - } else vj = mass[type[j]] / rho[j]; + rhoj = rho[j]; + if (interface_flag) + if (!(status[j] & STATUS_FLUID)) + rhoj = compute_interface->correct_rho(j,i); + vj = mass[type[j]] / rhoj; M += w * vj; } } @@ -596,11 +576,12 @@ void ComputeRHEOKernel::compute_peratom() r = sqrt(rsq); w = calc_w_quintic(i,j,dx[0],dx[1],dx[2],r); - if (solid_flag) + rhoj = rho[j]; + if (interface_flag) if (!(status[j] & STATUS_FLUID)) - vj = mass[type[j]]/compute_interface->correct_rho(j,i); - else - vj = mass[type[j]]/rho[j]; + rhoj = compute_interface->correct_rho(j,i); + + vj = mass[type[j]] / rhoj; //Populate the H-vector of polynomials (2D) if (dim == 2) { @@ -759,7 +740,7 @@ void ComputeRHEOKernel::compute_coordination() void ComputeRHEOKernel::grow_arrays(int nmax) { - memory->grow(coordination, nmax, "atom:rheo_coordination"); + memory->grow(coordination, nmax, "rheo:coordination"); if (kernel_style == CRK0) { memory->grow(C0, nmax, "rheo/kernel:C0"); diff --git a/src/RHEO/compute_rheo_kernel.h b/src/RHEO/compute_rheo_kernel.h index 1842406977..5324199f76 100644 --- a/src/RHEO/compute_rheo_kernel.h +++ b/src/RHEO/compute_rheo_kernel.h @@ -49,7 +49,7 @@ class ComputeRHEOKernel : public Compute { private: int comm_stage, comm_forward_save; - int solid_flag; + int interface_flag; int gsl_error_flag; std::unordered_set gsl_error_tags; diff --git a/src/RHEO/compute_rheo_property_atom.cpp b/src/RHEO/compute_rheo_property_atom.cpp index a23ba5a639..7682552abe 100644 --- a/src/RHEO/compute_rheo_property_atom.cpp +++ b/src/RHEO/compute_rheo_property_atom.cpp @@ -21,20 +21,29 @@ #include "atom.h" #include "atom_vec.h" +#include "compute_rheo_interface.h" +#include "compute_rheo_kernel.h" +#include "compute_rheo_surface.h" +#include "compute_rheo_vshift.h" #include "error.h" +#include "fix_rheo.h" +#include "fix_rheo_thermal.h" #include "memory.h" +#include "modify.h" #include "update.h" #include #include using namespace LAMMPS_NS; +using namespace RHEO_NS; /* ---------------------------------------------------------------------- */ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - index(nullptr), colindex(nullptr), pack_choice(nullptr) + Compute(lmp, narg, arg), fix_rheo(nullptr), fix_thermal(nullptr), compute_interface(nullptr), + compute_kernel(nullptr), compute_surface(nullptr), compute_vshift(nullptr), + index(nullptr), pack_choice(nullptr) { if (narg < 4) utils::missing_cmd_args(FLERR, "compute property/atom", error); @@ -43,31 +52,66 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a if (nvalues == 1) size_peratom_cols = 0; else size_peratom_cols = nvalues; + thermal_flag, interface_flag, surface_flag, shift_flag = 0; + // parse input values // customize a new keyword by adding to if statement pack_choice = new FnPtrPack[nvalues]; + index = new int[nvalues]; int i; for (int iarg = 3; iarg < narg; iarg++) { i = iarg-3; - if (strcmp(arg[iarg],"id") == 0) { - pack_choice[i] = &ComputeRHEOPropertyAtom::pack_id; - } else if (strcmp(arg[iarg],"mol") == 0) { - if (!atom->molecule_flag) - error->all(FLERR,"Compute property/atom {} is not available", arg[iarg]); - pack_choice[i] = &ComputeRHEOPropertyAtom::pack_molecule; - } else if (strcmp(arg[iarg],"proc") == 0) { - pack_choice[i] = &ComputeRHEOPropertyAtom::pack_proc; - } else if (strcmp(arg[iarg],"type") == 0) { - pack_choice[i] = &ComputeRHEOPropertyAtom::pack_type; - } else if (strcmp(arg[iarg],"mass") == 0) { - pack_choice[i] = &ComputeRHEOPropertyAtom::pack_mass; - - + if (strcmp(arg[iarg],"phase") == 0) { + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_phase; + } else if (strcmp(arg[iarg],"chi") == 0) { + interface_flag = 1; + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_chi; + } else if (strcmp(arg[iarg],"surface") == 0) { + surface_flag = 1; + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface; + } else if (strcmp(arg[iarg],"surface/r") == 0) { + surface_flag = 1; + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_r; + } else if (strcmp(arg[iarg],"surface/divr") == 0) { + surface_flag = 1; + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_divr; + } else if (strcmp(arg[iarg],"surface/nx") == 0) { + surface_flag = 1; + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_nx; + } else if (strcmp(arg[iarg],"surface/ny") == 0) { + surface_flag = 1; + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_ny; + } else if (strcmp(arg[iarg],"surface/nz") == 0) { + surface_flag = 1; + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_nz; + } else if (strcmp(arg[iarg],"coordination") == 0) { + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_coordination; + } else if (strcmp(arg[iarg],"cv") == 0) { + thermal_flag = 1; + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_cv; + } else if (strcmp(arg[iarg],"shift/vx") == 0) { + shift_flag = 1; + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_shift_vx; + } else if (strcmp(arg[iarg],"shift/vy") == 0) { + shift_flag = 1; + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_shift_vy; + } else if (strcmp(arg[iarg],"shift/vz") == 0) { + shift_flag = 1; + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_shift_vz; } else { - error->all(FLERR,"Invalid keyword {} for compute rheo/property/atom command ", arg[iarg]); + index[i] = atom->avec->property_atom(arg[iarg]); + if (index[i] < 0) + error->all(FLERR, + "Invalid keyword {} for atom style {} in compute rheo/property/atom command ", + atom->get_style(), arg[iarg]); + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_atom_style; + + if (strcmp(arg[iarg],"temperature") == 0) thermal_flag = 1; + if (strcmp(arg[iarg],"heatflow") == 0) thermal_flag = 1; + if (strcmp(arg[iarg],"conductivity") == 0) thermal_flag = 1; } } @@ -79,12 +123,41 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a ComputeRHEOPropertyAtom::~ComputeRHEOPropertyAtom() { delete[] pack_choice; + delete[] index; memory->destroy(vector_atom); memory->destroy(array_atom); } /* ---------------------------------------------------------------------- */ +void ComputeRHEOPropertyAtom::init() +{ + auto fixes = modify->get_fix_by_style("rheo"); + if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use fix rheo/viscosity"); + fix_rheo = dynamic_cast(fixes[0]); + + if (interface_flag && !(fix_rheo->interface_flag)) + error->all(FLERR, "Cannot request interfacial property without corresponding option in fix rheo"); + if (surface_flag && !(fix_rheo->surface_flag)) + error->all(FLERR, "Cannot request surface property without corresponding option in fix rheo"); + if (shift_flag && !(fix_rheo->shift_flag)) + error->all(FLERR, "Cannot request velocity shifting property without corresponding option in fix rheo"); + if (thermal_flag && !(fix_rheo->thermal_flag)) + error->all(FLERR, "Cannot request thermal property without fix rheo/thermal"); + + compute_interface = fix_rheo->compute_interface; + compute_kernel = fix_rheo->compute_kernel; + compute_surface = fix_rheo->compute_surface; + compute_vshift = fix_rheo->compute_vshift; + + if (thermal_flag) { + fixes = modify->get_fix_by_style("rheo/thermal"); + fix_thermal = dynamic_cast(fixes[0]); + } +} + +/* ---------------------------------------------------------------------- */ + void ComputeRHEOPropertyAtom::compute_peratom() { invoked_peratom = update->ntimestep; @@ -133,14 +206,16 @@ double ComputeRHEOPropertyAtom::memory_usage() /* ---------------------------------------------------------------------- */ -void ComputeRHEOPropertyAtom::pack_rho(int n) +void ComputeRHEOPropertyAtom::pack_phase(int n) { - double *rho = atom->rho; + int *status = atom->status; 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] = rho[i]; + if (mask[i] & groupbit) buf[n] = (status[i] & inverse_mask); else buf[n] = 0.0; n += nvalues; } @@ -148,15 +223,189 @@ void ComputeRHEOPropertyAtom::pack_rho(int n) /* ---------------------------------------------------------------------- */ -void ComputeRHEOPropertyAtom::pack_status(int n) +void ComputeRHEOPropertyAtom::pack_chi(int n) +{ + double *chi = compute_interface->chi; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = chi[i]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_surface(int n) { int *status = atom->status; 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]; + if (mask[i] & groupbit) buf[n] = (status[i] & inverse_mask); else buf[n] = 0.0; n += nvalues; } } + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_surface_r(int n) +{ + double *rsurface = compute_surface->rsurface; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = rsurface[i]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_surface_divr(int n) +{ + double *divr = compute_surface->divr; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = divr[i]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_surface_nx(int n) +{ + double **nsurface = compute_surface->nsurface; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = nsurface[i][0]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_surface_ny(int n) +{ + double **nsurface = compute_surface->nsurface; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = nsurface[i][1]; + else buf[n] = 0.0; + n += nvalues; + } +} + + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_surface_nz(int n) +{ + double **nsurface = compute_surface->nsurface; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = nsurface[i][2]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_coordination(int n) +{ + int *coordination = compute_kernel->coordination; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = coordination[i]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_cv(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = fix_thermal->calc_cv(i); + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_shift_vx(int n) +{ + double **vshift = compute_vshift->vshift; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = vshift[i][0]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_shift_vy(int n) +{ + double **vshift = compute_vshift->vshift; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = vshift[i][1]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_shift_vz(int n) +{ + double **vshift = compute_vshift->vshift; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = vshift[i][2]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_atom_style(int n) +{ + atom->avec->pack_property_atom(index[n],&buf[n],nvalues,groupbit); +} diff --git a/src/RHEO/compute_rheo_property_atom.h b/src/RHEO/compute_rheo_property_atom.h index 26ca004da0..accb7e9156 100644 --- a/src/RHEO/compute_rheo_property_atom.h +++ b/src/RHEO/compute_rheo_property_atom.h @@ -28,37 +28,40 @@ class ComputeRHEOPropertyAtom : public Compute { public: ComputeRHEOPropertyAtom(class LAMMPS *, int, char **); ~ComputeRHEOPropertyAtom() override; + void init() override; void compute_peratom() override; double memory_usage() override; private: - int nvalues; - int nmax; + int nvalues, nmax; + int thermal_flag, interface_flag, surface_flag, shift_flag; + int *index; double *buf; typedef void (ComputeRHEOPropertyAtom::*FnPtrPack)(int); FnPtrPack *pack_choice; // ptrs to pack functions - void pack_rho(int); - void pack_drho(int); - void pack_temperature(int); - void pack_heatflow(int); - void pack_status(int); void pack_phase(int); + void pack_chi(int); void pack_surface(int); - void pack_r_surface(int); - void pack_divr_surface(int); - void pack_nx_surface(int); - void pack_ny_surface(int); - void pack_nz_surface(int); + void pack_surface_r(int); + void pack_surface_divr(int); + void pack_surface_nx(int); + void pack_surface_ny(int); + void pack_surface_nz(int); void pack_coordination(int); - void pack_viscosity(int); - void pack_pressure(int); - void pack_conductivity(int); void pack_cv(int); - void pack_vx_shift(int); - void pack_vy_shift(int); - void pack_vz_shift(int); + void pack_shift_vx(int); + void pack_shift_vy(int); + void pack_shift_vz(int); + void pack_atom_style(int); + + class FixRHEO *fix_rheo; + class FixRHEOThermal *fix_thermal; + class ComputeRHEOInterface *compute_interface; + class ComputeRHEOKernel *compute_kernel; + class ComputeRHEOSurface *compute_surface; + class ComputeRHEOVShift *compute_vshift; }; diff --git a/src/RHEO/compute_rheo_surface.h b/src/RHEO/compute_rheo_surface.h index 58a3e3b9c4..220f8beb6d 100644 --- a/src/RHEO/compute_rheo_surface.h +++ b/src/RHEO/compute_rheo_surface.h @@ -36,13 +36,13 @@ class ComputeRHEOSurface : public Compute { int pack_forward_comm(int, int *, double *, int, int *) override; void unpack_forward_comm(int, int, double *) override; - double **nsurface, *rsurface; + double **nsurface, *rsurface, *divr; class FixRHEO *fix_rheo; private: double cut, cutsq, rho0, threshold_divr; int surface_style, nmax_store, threshold_z; - double **B, **gradC, *divr; + double **B, **gradC; int threshold_style, comm_stage; class NeighList *list; diff --git a/src/RHEO/compute_rheo_vshift.cpp b/src/RHEO/compute_rheo_vshift.cpp index 440bfc7fc7..3d3914436e 100644 --- a/src/RHEO/compute_rheo_vshift.cpp +++ b/src/RHEO/compute_rheo_vshift.cpp @@ -22,6 +22,7 @@ #include "comm.h" #include "compute_rheo_interface.h" #include "compute_rheo_kernel.h" +#include "compute_rheo_surface.h" #include "domain.h" #include "error.h" #include "fix_rheo.h" @@ -38,36 +39,22 @@ using namespace RHEO_NS; ComputeRHEOVShift::ComputeRHEOVShift(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), list(nullptr), vshift(nullptr), fix_rheo(nullptr), - compute_kernel(nullptr), compute_interface(nullptr) + compute_kernel(nullptr), compute_interface(nullptr), compute_surface(nullptr) { if (narg != 3) error->all(FLERR,"Illegal compute RHEO/VShift command"); comm_reverse = 3; surface_flag = 0; - // 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_store exceeded - - int tmp1, tmp2; - int index = atom->find_custom("rheo_vshift", tmp1, tmp2); - if (index == -1) { - index = atom->add_custom("rheo_vshift", 1, 3); - nmax_store = atom->nmax; - } - vshift = atom->darray[index]; + nmax_store = atom->nmax; + memory->create(vshift, nmax_store, 3, "rheo:vshift"); } /* ---------------------------------------------------------------------- */ ComputeRHEOVShift::~ComputeRHEOVShift() { - // 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, 1, 3); - + memory->destroy(vshift); } /* ---------------------------------------------------------------------- */ @@ -76,12 +63,12 @@ void ComputeRHEOVShift::init() { neighbor->add_request(this, NeighConst::REQ_DEFAULT); - surface_flag = 0; - if (fix_rheo->surface_flag) - surface_flag = 1; + surface_flag = fix_rheo->surface_flag; + interface_flag = fix_rheo->interface_flag; compute_kernel = fix_rheo->compute_kernel; compute_interface = fix_rheo->compute_interface; + compute_surface = fix_rheo->compute_surface; cut = fix_rheo->cut; cutsq = cut * cut; @@ -128,7 +115,7 @@ void ComputeRHEOVShift::compute_peratom() firstneigh = list->firstneigh; if (nmax_store < atom->nmax) { - memory->grow(vshift, atom->nmax, 3, "atom:rheo_vshift"); + memory->grow(vshift, atom->nmax, 3, "rheo:vshift"); nmax_store = atom->nmax; } @@ -176,15 +163,17 @@ void ComputeRHEOVShift::compute_peratom() rhoj = rho[j]; // Add corrections for walls - if (fluidi && (!fluidj)) { - compute_interface->correct_v(vi, vj, i, j); - rhoj = compute_interface->correct_rho(j,i); - } else if ((!fluidi) && fluidj) { - compute_interface->correct_v(vj, vi, j, i); - rhoi = compute_interface->correct_rho(i,j); - } else if ((!fluidi) && (!fluidj)) { - rhoi = 1.0; - rhoj = 1.0; + if (interface_flag) { + if (fluidi && (!fluidj)) { + compute_interface->correct_v(vi, vj, i, j); + rhoj = compute_interface->correct_rho(j,i); + } else if ((!fluidi) && fluidj) { + compute_interface->correct_v(vj, vi, j, i); + rhoi = compute_interface->correct_rho(i,j); + } else if ((!fluidi) && (!fluidj)) { + rhoi = 1.0; + rhoj = 1.0; + } } voli = imass / rhoi; @@ -227,30 +216,28 @@ void ComputeRHEOVShift::correct_surfaces() { if (!surface_flag) return; + int i, a, b; + int *status = atom->status; int *mask = atom->mask; - int nlocal = atom->nlocal; - int i, a, b; - int dim = domain->dimension; + double **nsurface = compute_surface->nsurface; - 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]; + int nlocal = atom->nlocal; + int dim = domain->dimension; double nx,ny,nz,vx,vy,vz; for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { if ((status[i] & STATUS_SURFACE) || (status[i] & STATUS_LAYER)) { - nx = nsurf[i][0]; - ny = nsurf[i][1]; + nx = nsurface[i][0]; + ny = nsurface[i][1]; vx = vshift[i][0]; vy = vshift[i][1]; vz = vshift[i][2]; vshift[i][0] = (1 - nx * nx) * vx - nx * ny * vy; vshift[i][1] = (1 - ny * ny) * vy - nx * ny * vx; if (dim > 2) { - nz = nsurf[i][2]; + nz = nsurface[i][2]; vshift[i][0] -= nx * nz * vz; vshift[i][1] -= ny * nz * vz; vshift[i][2] = (1 - nz * nz) * vz - nz * ny * vy - nx * nz * vx; diff --git a/src/RHEO/compute_rheo_vshift.h b/src/RHEO/compute_rheo_vshift.h index 8611e177d1..9d3a0166d6 100644 --- a/src/RHEO/compute_rheo_vshift.h +++ b/src/RHEO/compute_rheo_vshift.h @@ -42,11 +42,12 @@ class ComputeRHEOVShift : public Compute { private: int nmax_store; double dtv, cut, cutsq, cutthird; - int surface_flag; + int surface_flag, interface_flag; class NeighList *list; - class ComputeRHEOInterface *compute_interface ; + class ComputeRHEOInterface *compute_interface; class ComputeRHEOKernel *compute_kernel; + class ComputeRHEOSurface *compute_surface; }; } // namespace LAMMPS_NS diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index ff804fe007..ac870affd5 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -60,8 +60,12 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : if (igroup != 0) error->all(FLERR,"fix rheo command requires group all"); + if (atom->pressure_flag != 1) + error->all(FLERR,"fix rheo command requires atom_style with pressure"); if (atom->rho_flag != 1) error->all(FLERR,"fix rheo command requires atom_style with density"); + if (atom->viscosity_flag != 1) + error->all(FLERR,"fix rheo command requires atom_style with viscosity"); if (atom->status_flag != 1) error->all(FLERR,"fix rheo command requires atom_style with status"); diff --git a/src/RHEO/fix_rheo_pressure.cpp b/src/RHEO/fix_rheo_pressure.cpp index 75edf42572..45557794cd 100644 --- a/src/RHEO/fix_rheo_pressure.cpp +++ b/src/RHEO/fix_rheo_pressure.cpp @@ -38,7 +38,7 @@ static constexpr double SEVENTH = 1.0 / 7.0; /* ---------------------------------------------------------------------- */ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), fix_rheo(nullptr), pressure(nullptr) + Fix(lmp, narg, arg), fix_rheo(nullptr) { if (narg < 4) error->all(FLERR,"Illegal fix command"); diff --git a/src/RHEO/fix_rheo_thermal.cpp b/src/RHEO/fix_rheo_thermal.cpp index 5df8c1c506..0d8909b908 100644 --- a/src/RHEO/fix_rheo_thermal.cpp +++ b/src/RHEO/fix_rheo_thermal.cpp @@ -40,7 +40,7 @@ enum {NONE, CONSTANT, TYPE}; FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), fix_rheo(nullptr), compute_grad(nullptr), compute_vshift(nullptr), - Tc_type(nullptr), kappa_type(nullptr), cv_type(nullptr), conductivity(nullptr) + Tc_type(nullptr), kappa_type(nullptr), cv_type(nullptr) { if (narg < 4) error->all(FLERR,"Illegal fix command"); @@ -48,9 +48,6 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) : cv_style = NONE; conductivity_style = NONE; - comm_forward = 1; - nmax_store = 0; - int ntypes = atom->ntypes; int iarg = 3; while (iarg < narg) { @@ -170,9 +167,11 @@ void FixRHEOThermal::init() dtf = 0.5 * update->dt * force->ftm2v; if (atom->temperature_flag != 1) - error->all(FLERR,"fix rheo/thermal command requires atoms store temperature property"); + error->all(FLERR,"fix rheo/thermal command requires atom property temperature"); if (atom->heatflow_flag != 1) - error->all(FLERR,"fix rheo/thermal command requires atoms store heatflow property"); + error->all(FLERR,"fix rheo/thermal command requires atom property heatflow"); + if (atom->conductivity_flag != 1) + error->all(FLERR,"fix rheo/thermal command requires atom property conductivity"); } /* ---------------------------------------------------------------------- */ @@ -181,34 +180,6 @@ void FixRHEOThermal::setup_pre_force(int /*vflag*/) { fix_rheo->thermal_fix_defined = 1; - // Identify whether this is the first/last instance of fix thermal - // First will grow arrays, last will communicate - first_flag = 0; - last_flag = 0; - - int i = 0; - auto fixlist = modify->get_fix_by_style("rheo/thermal"); - for (const auto &fix : fixlist) { - if (strcmp(fix->id, id) == 0) break; - i++; - } - - if (i == 0) first_flag = 1; - if ((i + 1) == fixlist.size()) last_flag = 1; - - // Create conductivity 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_conductivity", tmp1, tmp2); - if (index== -1) { - index = atom->add_custom("rheo_conductivity", 1, 0); - nmax_store = atom->nmax; - } - conductivity = atom->dvector[index]; - post_neighbor(); pre_force(0); } @@ -294,14 +265,10 @@ void FixRHEOThermal::post_neighbor() int i; int *type = atom->type; int *mask = atom->mask; + double *conductivity = atom->conductivity; int nlocal = atom->nlocal; int nall = nlocal + atom->nghost; - if (first_flag && (nmax_store < atom->nmax)) { - memory->grow(conductivity, atom->nmax, "atom:rheo_conductivity"); - nmax_store = atom->nmax; - } - if (conductivity_style == CONSTANT) { for (i = 0; i < nall; i++) if (mask[i] & groupbit) conductivity[i] = kappa; @@ -312,39 +279,11 @@ void FixRHEOThermal::post_neighbor() } /* ---------------------------------------------------------------------- - Update (and forward) evolving conductivity styles every timestep - Zero heat flow + In the future, update & forward evolving conductivity styles every timestep ------------------------------------------------------------------------- */ void FixRHEOThermal::pre_force(int /*vflag*/) { - // send updated temperatures to ghosts if first instance of fix - // then clear heatflow for next force calculation - double *heatflow = atom->heatflow; - if (first_flag) { - comm->forward_comm(this); - for (int i = 0; i < atom->nmax; i++) heatflow[i] = 0.0; - } - - // Not needed yet, when needed add stage check for (un)pack_forward_comm() methods - //int i; - //double *conductivity = atom->dvector[index_cond]; - //int *mask = atom->mask; - //int nlocal = atom->nlocal; - - //if (first_flag && (nmax_store < atom->nmax)) { - // memory->grow(conductivity, atom->nmax, "atom:rheo_conductivity"); - // nmax_store = atom->nmax; - //} - - //if (conductivity_style == TBD) { - // for (i = 0; i < nlocal; i++) { - // if (mask[i] & groupbit) { - // } - // } - //} - - //if (last_flag && comm_forward) comm->forward_comm(this); } /* ---------------------------------------------------------------------- */ @@ -387,69 +326,3 @@ double FixRHEOThermal::calc_cv(int i) return(cv_type[atom->type[i]]); } } - -/* ---------------------------------------------------------------------- */ - -int FixRHEOThermal::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/) -{ - int i, j, m; - - double *temperature = atom->temperature; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = temperature[j]; - } - - return m; -} - -/* ---------------------------------------------------------------------- */ - -void FixRHEOThermal::unpack_forward_comm(int n, int first, double *buf) -{ - int i, m, last; - - m = 0; - last = first + n; - - double *temperature = atom->temperature; - - for (i = first; i < last; i++) temperature[i] = buf[m++]; -} - -/* ---------------------------------------------------------------------- */ - -int FixRHEOThermal::pack_reverse_comm(int n, int first, double *buf) -{ - int m = 0; - int last = first + n; - double *heatflow = atom->heatflow; - - for (int i = first; i < last; i++) { - buf[m++] = heatflow[i]; - } - - return m; -} - -/* ---------------------------------------------------------------------- */ - -void FixRHEOThermal::unpack_reverse_comm(int n, int *list, double *buf) -{ - int m = 0; - double *heatflow = atom->heatflow; - - for (int i = 0; i < n; i++) - heatflow[list[i]] += buf[m++]; -} - -/* ---------------------------------------------------------------------- */ - -double FixRHEOThermal::memory_usage() -{ - double bytes = 0.0; - bytes += (size_t) nmax_store * sizeof(double); - return bytes; -} diff --git a/src/RHEO/fix_rheo_thermal.h b/src/RHEO/fix_rheo_thermal.h index cf64c0b8d1..a27ad98a8c 100644 --- a/src/RHEO/fix_rheo_thermal.h +++ b/src/RHEO/fix_rheo_thermal.h @@ -37,29 +37,22 @@ class FixRHEOThermal : public Fix { void pre_force(int) override; void final_integrate() override; void reset_dt() 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; + double calc_cv(int); private: double *cv_type, cv; double *Tc_type, Tc; double *kappa_type, kappa; double dtf, dtv; - double *conductivity; int Tc_style; int cv_style; int conductivity_style; - int first_flag, last_flag; - int nmax_store; class FixRHEO *fix_rheo; class ComputeRHEOGrad *compute_grad; class ComputeRHEOVShift *compute_vshift; - double calc_cv(int); + void grow_array(int); }; } // namespace LAMMPS_NS diff --git a/src/RHEO/fix_rheo_viscosity.cpp b/src/RHEO/fix_rheo_viscosity.cpp index b15f488370..4d70ffaca3 100644 --- a/src/RHEO/fix_rheo_viscosity.cpp +++ b/src/RHEO/fix_rheo_viscosity.cpp @@ -37,14 +37,13 @@ enum {NONE, CONSTANT, TYPE, POWER}; /* ---------------------------------------------------------------------- */ FixRHEOViscosity::FixRHEOViscosity(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), fix_rheo(nullptr), eta_type(nullptr), viscosity(nullptr), compute_grad(nullptr) + Fix(lmp, narg, arg), fix_rheo(nullptr), compute_grad(nullptr), eta_type(nullptr) { if (narg < 4) error->all(FLERR,"Illegal fix command"); viscosity_style = NONE; comm_forward = 0; - nmax_store = 0; int ntypes = atom->ntypes; int iarg = 3; @@ -86,11 +85,6 @@ FixRHEOViscosity::FixRHEOViscosity(LAMMPS *lmp, int narg, char **arg) : FixRHEOViscosity::~FixRHEOViscosity() { - // Remove custom property if it exists - int tmp1, tmp2, index; - index = atom->find_custom("rheo_viscosity", tmp1, tmp2); - if (index != -1) atom->remove_custom(index, 1, 0); - memory->destroy(eta_type); } @@ -121,9 +115,7 @@ void FixRHEOViscosity::setup_pre_force(int /*vflag*/) { fix_rheo->viscosity_fix_defined = 1; - // Identify whether this is the first/last instance of fix viscosity - // First will grow arrays, last will communicate - first_flag = 0; + // Identify whether this is the last instance of fix viscosity last_flag = 0; int i = 0; @@ -133,22 +125,8 @@ void FixRHEOViscosity::setup_pre_force(int /*vflag*/) i++; } - if (i == 0) first_flag = 1; if ((i + 1) == fixlist.size()) last_flag = 1; - // Create viscosity 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_viscosity", tmp1, tmp2); - if (index == -1) { - index = atom->add_custom("rheo_viscosity", 1, 0); - nmax_store = atom->nmax; - } - viscosity = atom->dvector[index]; - post_neighbor(); pre_force(0); } @@ -163,14 +141,9 @@ void FixRHEOViscosity::post_neighbor() int *type = atom->type; int *mask = atom->mask; + double *viscosity = atom->viscosity; - int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; - - if (first_flag && (nmax_store < atom->nmax)) { - memory->grow(viscosity, atom->nmax, "atom:rheo_viscosity"); - nmax_store = atom->nmax; - } + int nall = atom->nlocal + atom->nghost; if (viscosity_style == CONSTANT) { for (i = 0; i < nall; i++) @@ -191,16 +164,12 @@ void FixRHEOViscosity::pre_force(int /*vflag*/) double tmp, gdot; int *mask = atom->mask; + double *viscosity = atom->viscosity; double **gradv = compute_grad->gradv; int nlocal = atom->nlocal; int dim = domain->dimension; - if (first_flag && (nmax_store < atom->nmax)) { - memory->grow(viscosity, atom->nmax, "atom:rheo_viscosity"); - nmax_store = atom->nmax; - } - if (viscosity_style == POWER) { for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { @@ -231,7 +200,8 @@ void FixRHEOViscosity::pre_force(int /*vflag*/) int FixRHEOViscosity::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/) { - int i,j,k,m; + int i, j, k, m; + double *viscosity = atom->viscosity; m = 0; for (i = 0; i < n; i++) { @@ -246,6 +216,7 @@ int FixRHEOViscosity::pack_forward_comm(int n, int *list, double *buf, void FixRHEOViscosity::unpack_forward_comm(int n, int first, double *buf) { int i, k, m, last; + double *viscosity = atom->viscosity; m = 0; last = first + n; @@ -253,12 +224,3 @@ void FixRHEOViscosity::unpack_forward_comm(int n, int first, double *buf) viscosity[i] = buf[m++]; } } - -/* ---------------------------------------------------------------------- */ - -double FixRHEOViscosity::memory_usage() -{ - double bytes = 0.0; - bytes += (size_t) nmax_store * sizeof(double); - return bytes; -} diff --git a/src/RHEO/fix_rheo_viscosity.h b/src/RHEO/fix_rheo_viscosity.h index 66df51601e..c681d18c00 100644 --- a/src/RHEO/fix_rheo_viscosity.h +++ b/src/RHEO/fix_rheo_viscosity.h @@ -35,15 +35,11 @@ class FixRHEOViscosity : 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; private: double *eta_type, eta; double npow, K, gd0, tau0; - double *viscosity; - int viscosity_style; - int first_flag, last_flag; - int nmax_store; + int viscosity_style, last_flag; class FixRHEO *fix_rheo; class ComputeRHEOGrad *compute_grad; diff --git a/src/RHEO/pair_rheo.cpp b/src/RHEO/pair_rheo.cpp index c61d613d82..03e859316f 100644 --- a/src/RHEO/pair_rheo.cpp +++ b/src/RHEO/pair_rheo.cpp @@ -99,6 +99,9 @@ void PairRHEO::compute(int eflag, int vflag) double *rho = atom->rho; double *mass = atom->mass; double *drho = atom->drho; + double *pressure = atom->pressure; + double *viscosity = atom->viscosity; + double *conductivity = atom->conductivity; double *temperature = atom->temperature; double *heatflow = atom->heatflow; double *special_lj = force->special_lj; @@ -112,22 +115,6 @@ void PairRHEO::compute(int eflag, int vflag) chi = compute_interface->chi; } - int tmp1, tmp2; - int index = atom->find_custom("rheo_viscosity", tmp1, tmp2); - if (index == -1) error->all(FLERR, "Cannot find rheo viscosity"); - double *viscosity = atom->dvector[index]; - - index = atom->find_custom("rheo_pressure", tmp1, tmp2); - if (index == -1) error->all(FLERR, "Cannot find rheo pressure"); - double *pressure = atom->dvector[index]; - - double *conductivity; - if (thermal_flag) { - index = atom->find_custom("rheo_conductivity", tmp1, tmp2); - if (index == -1) error->all(FLERR, "Cannot find rheo conductivity"); - conductivity = atom->dvector[index]; - } - inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; diff --git a/src/atom.cpp b/src/atom.cpp index 25023e0d49..cf8cb7468a 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -200,6 +200,9 @@ Atom::Atom(LAMMPS *_lmp) : Pointers(_lmp) // RHEO package status = nullptr; + conductivity = nullptr; + pressure = nullptr; + viscosity = nullptr; // SPH package @@ -533,6 +536,9 @@ void Atom::peratom_create() // RHEO package add_peratom("status",&status,INT,0); + add_peratom("conductivity",&conductivity,DOUBLE,0); + add_peratom("pressure",&pressure,DOUBLE,0); + add_peratom("viscosity",&viscosity,DOUBLE,0); // SPH package @@ -639,7 +645,7 @@ void Atom::set_atomflag_defaults() temperature_flag = heatflow_flag = 0; vfrac_flag = spin_flag = eradius_flag = ervel_flag = erforce_flag = 0; cs_flag = csforce_flag = vforce_flag = ervelforce_flag = etag_flag = 0; - status_flag = 0; + status_flag = conductivity_flag = pressure_flag = viscosity_flag = 0; rho_flag = esph_flag = cv_flag = vest_flag = 0; dpd_flag = edpd_flag = tdpd_flag = 0; sp_flag = 0; @@ -2943,6 +2949,9 @@ void *Atom::extract(const char *name) // RHEO package if (strcmp(name,"status") == 0) return (void *) status; + if (strcmp(name,"conductivity") == 0) return (void *) conductivity; + if (strcmp(name,"pressure") == 0) return (void *) pressure; + if (strcmp(name,"viscosity") == 0) return (void *) viscosity; // SPH package @@ -3069,6 +3078,9 @@ int Atom::extract_datatype(const char *name) // RHEO package if (strcmp(name,"status") == 0) return LAMMPS_INT; + if (strcmp(name,"conductivity") == 0) return LAMMPS_DOUBLE; + if (strcmp(name,"pressure") == 0) return LAMMPS_DOUBLE; + if (strcmp(name,"viscosity") == 0) return LAMMPS_DOUBLE; // SPH package diff --git a/src/atom.h b/src/atom.h index 51665cb0bf..a8ced43b4d 100644 --- a/src/atom.h +++ b/src/atom.h @@ -158,6 +158,9 @@ class Atom : protected Pointers { // RHEO package int *status; + double *conductivity; + double *pressure; + double *viscosity; // SPH package @@ -194,7 +197,7 @@ class Atom : protected Pointers { int temperature_flag, heatflow_flag; int vfrac_flag, spin_flag, eradius_flag, ervel_flag, erforce_flag; int cs_flag, csforce_flag, vforce_flag, ervelforce_flag, etag_flag; - int status_flag; + int status_flag, conductivity_flag, pressure_flag, viscosity_flag; int rho_flag, esph_flag, cv_flag, vest_flag; int dpd_flag, edpd_flag, tdpd_flag; int mesont_flag; diff --git a/src/atom_vec_rheo_thermal.cpp b/src/atom_vec_rheo_thermal.cpp new file mode 100644 index 0000000000..de0c7fa5d7 --- /dev/null +++ b/src/atom_vec_rheo_thermal.cpp @@ -0,0 +1,200 @@ +/* ---------------------------------------------------------------------- + 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 "atom_vec_rheo_thermal.h" + +#include "atom.h" + +#include + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +AtomVecRHEOThermal::AtomVecRHEOThermal(LAMMPS *lmp) : AtomVec(lmp) +{ + molecular = Atom::ATOMIC; + mass_type = PER_TYPE; + forceclearflag = 1; + + atom->status_flag = 1; + atom->conductivity_flag = 1; + atom->temperature_flag = 1; + atom->heatflow_flag = 1; + atom->pressure_flag = 1; + atom->rho_flag = 1; + atom->viscosity_flag = 1; + + // strings with peratom variables to include in each AtomVec method + // strings cannot contain fields in corresponding AtomVec default strings + // order of fields in a string does not matter + // except: fields_data_atom & fields_data_vel must match data file + + fields_grow = {"status", "rho", "drho", "temperature", "heatflow", "conductivity", "pressure", "viscosity"}; + fields_copy = {"status", "rho", "drho", "temperature", "heatflow", "conductivity", "pressure", "viscosity"}; + fields_comm = {"status", "rho", "temperature"}; + fields_comm_vel = {"status", "rho", "temperature"}; + fields_reverse = {"drho", "heatflow"}; + fields_border = {"status", "rho", "temperature"}; + fields_border_vel = {"status", "rho", "temperature"}; + fields_exchange = {"status", "rho", "temperature"}; + fields_restart = {"status", "rho", "temperature"}; + fields_create = {"status", "rho", "drho", "temperature", "heatflow", "conductivity", "pressure", "viscosity"}; + fields_data_atom = {"id", "type", "status", "rho", "temperature", "x"}; + fields_data_vel = {"id", "v"}; + + setup_fields(); +} + +/* ---------------------------------------------------------------------- + set local copies of all grow ptrs used by this class, except defaults + needed in replicate when 2 atom classes exist and it calls pack_restart() +------------------------------------------------------------------------- */ + +void AtomVecRHEOThermal::grow_pointers() +{ + status = atom->status; + conductivity = atom->conductivity; + temperature = atom->temperature; + heatflow = atom->heatflow; + pressure = atom->pressure; + rho = atom->rho; + drho = atom->drho; + viscosity = atom->viscosity; +} + +/* ---------------------------------------------------------------------- + clear extra forces starting at atom N + nbytes = # of bytes to clear for a per-atom vector +------------------------------------------------------------------------- */ + +void AtomVecRHEOThermal::force_clear(int n, size_t nbytes) +{ + memset(&drho[n], 0, nbytes); + memset(&heatflow[n], 0, nbytes); +} + +/* ---------------------------------------------------------------------- + modify what AtomVec::data_atom() just unpacked + or initialize other atom quantities +------------------------------------------------------------------------- */ + +void AtomVecRHEOThermal::data_atom_post(int ilocal) +{ + drho[ilocal] = 0.0; + heatflow[ilocal] = 0.0; + pressure[ilocal] = 0.0; + viscosity[ilocal] = 0.0; + conductivity[ilocal] = 0.0; +} + +/* ---------------------------------------------------------------------- + assign an index to named atom property and return index + return -1 if name is unknown to this atom style +------------------------------------------------------------------------- */ + +int AtomVecRHEOThermal::property_atom(const std::string &name) +{ + if (name == "status") return 0; + if (name == "rho") return 1; + if (name == "drho") return 2; + if (name == "temperature") return 3; + if (name == "heatflow") return 4; + if (name == "conductivity") return 5; + if (name == "pressure") return 6; + if (name == "viscosity") return 7; + return -1; +} + +/* ---------------------------------------------------------------------- + pack per-atom data into buf for ComputePropertyAtom + index maps to data specific to this atom style +------------------------------------------------------------------------- */ + +void AtomVecRHEOThermal::pack_property_atom(int index, double *buf, int nvalues, int groupbit) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + int n = 0; + + if (index == 0) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = status[i]; + else + buf[n] = 0.0; + n += nvalues; + } + } else if (index == 1) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = rho[i]; + else + buf[n] = 0.0; + n += nvalues; + } + } else if (index == 2) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = drho[i]; + else + buf[n] = 0.0; + n += nvalues; + } + } else if (index == 3) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = temperature[i]; + else + buf[n] = 0.0; + n += nvalues; + } + } else if (index == 4) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = heatflow[i]; + else + buf[n] = 0.0; + n += nvalues; + } + } else if (index == 5) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = conductivity[i]; + else + buf[n] = 0.0; + n += nvalues; + } + } else if (index == 6) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = pressure[i]; + else + buf[n] = 0.0; + n += nvalues; + } + } else if (index == 7) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = viscosity[i]; + else + buf[n] = 0.0; + n += nvalues; + } + } +} diff --git a/src/atom_vec_rheo_thermal.h b/src/atom_vec_rheo_thermal.h new file mode 100644 index 0000000000..27c6c3c9b5 --- /dev/null +++ b/src/atom_vec_rheo_thermal.h @@ -0,0 +1,46 @@ +/* -*- 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 ATOM_CLASS +// clang-format off +AtomStyle(rheo/thermal,AtomVecRHEOThermal); +// clang-format on +#else + +#ifndef LMP_ATOM_VEC_RHEO_THERMAL_H +#define LMP_ATOM_VEC_RHEO_THERMAL_H + +#include "atom_vec.h" + +namespace LAMMPS_NS { + +class AtomVecRHEOThermal : virtual public AtomVec { + public: + AtomVecRHEOThermal(class LAMMPS *); + + void grow_pointers() override; + void force_clear(int, size_t) override; + void data_atom_post(int) override; + int property_atom(const std::string &) override; + void pack_property_atom(int, double *, int, int) override; + + private: + int *status; + double *conductivity, *temperature, *heatflow; + double *pressure, *rho, *drho, *viscosity; +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/compute_rheo_property_atom.cpp b/src/compute_rheo_property_atom.cpp new file mode 100644 index 0000000000..7682552abe --- /dev/null +++ b/src/compute_rheo_property_atom.cpp @@ -0,0 +1,411 @@ +// clang-format off +/* ---------------------------------------------------------------------- + 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_property_atom.h" + +#include "atom.h" +#include "atom_vec.h" +#include "compute_rheo_interface.h" +#include "compute_rheo_kernel.h" +#include "compute_rheo_surface.h" +#include "compute_rheo_vshift.h" +#include "error.h" +#include "fix_rheo.h" +#include "fix_rheo_thermal.h" +#include "memory.h" +#include "modify.h" +#include "update.h" + +#include +#include + +using namespace LAMMPS_NS; +using namespace RHEO_NS; + +/* ---------------------------------------------------------------------- */ + +ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), fix_rheo(nullptr), fix_thermal(nullptr), compute_interface(nullptr), + compute_kernel(nullptr), compute_surface(nullptr), compute_vshift(nullptr), + index(nullptr), pack_choice(nullptr) +{ + if (narg < 4) utils::missing_cmd_args(FLERR, "compute property/atom", error); + + peratom_flag = 1; + nvalues = narg - 3; + if (nvalues == 1) size_peratom_cols = 0; + else size_peratom_cols = nvalues; + + thermal_flag, interface_flag, surface_flag, shift_flag = 0; + + // parse input values + // customize a new keyword by adding to if statement + + pack_choice = new FnPtrPack[nvalues]; + index = new int[nvalues]; + + int i; + for (int iarg = 3; iarg < narg; iarg++) { + i = iarg-3; + + if (strcmp(arg[iarg],"phase") == 0) { + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_phase; + } else if (strcmp(arg[iarg],"chi") == 0) { + interface_flag = 1; + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_chi; + } else if (strcmp(arg[iarg],"surface") == 0) { + surface_flag = 1; + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface; + } else if (strcmp(arg[iarg],"surface/r") == 0) { + surface_flag = 1; + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_r; + } else if (strcmp(arg[iarg],"surface/divr") == 0) { + surface_flag = 1; + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_divr; + } else if (strcmp(arg[iarg],"surface/nx") == 0) { + surface_flag = 1; + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_nx; + } else if (strcmp(arg[iarg],"surface/ny") == 0) { + surface_flag = 1; + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_ny; + } else if (strcmp(arg[iarg],"surface/nz") == 0) { + surface_flag = 1; + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_nz; + } else if (strcmp(arg[iarg],"coordination") == 0) { + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_coordination; + } else if (strcmp(arg[iarg],"cv") == 0) { + thermal_flag = 1; + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_cv; + } else if (strcmp(arg[iarg],"shift/vx") == 0) { + shift_flag = 1; + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_shift_vx; + } else if (strcmp(arg[iarg],"shift/vy") == 0) { + shift_flag = 1; + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_shift_vy; + } else if (strcmp(arg[iarg],"shift/vz") == 0) { + shift_flag = 1; + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_shift_vz; + } else { + index[i] = atom->avec->property_atom(arg[iarg]); + if (index[i] < 0) + error->all(FLERR, + "Invalid keyword {} for atom style {} in compute rheo/property/atom command ", + atom->get_style(), arg[iarg]); + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_atom_style; + + if (strcmp(arg[iarg],"temperature") == 0) thermal_flag = 1; + if (strcmp(arg[iarg],"heatflow") == 0) thermal_flag = 1; + if (strcmp(arg[iarg],"conductivity") == 0) thermal_flag = 1; + } + } + + nmax = 0; +} + +/* ---------------------------------------------------------------------- */ + +ComputeRHEOPropertyAtom::~ComputeRHEOPropertyAtom() +{ + delete[] pack_choice; + delete[] index; + memory->destroy(vector_atom); + memory->destroy(array_atom); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::init() +{ + auto fixes = modify->get_fix_by_style("rheo"); + if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use fix rheo/viscosity"); + fix_rheo = dynamic_cast(fixes[0]); + + if (interface_flag && !(fix_rheo->interface_flag)) + error->all(FLERR, "Cannot request interfacial property without corresponding option in fix rheo"); + if (surface_flag && !(fix_rheo->surface_flag)) + error->all(FLERR, "Cannot request surface property without corresponding option in fix rheo"); + if (shift_flag && !(fix_rheo->shift_flag)) + error->all(FLERR, "Cannot request velocity shifting property without corresponding option in fix rheo"); + if (thermal_flag && !(fix_rheo->thermal_flag)) + error->all(FLERR, "Cannot request thermal property without fix rheo/thermal"); + + compute_interface = fix_rheo->compute_interface; + compute_kernel = fix_rheo->compute_kernel; + compute_surface = fix_rheo->compute_surface; + compute_vshift = fix_rheo->compute_vshift; + + if (thermal_flag) { + fixes = modify->get_fix_by_style("rheo/thermal"); + fix_thermal = dynamic_cast(fixes[0]); + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::compute_peratom() +{ + invoked_peratom = update->ntimestep; + + // grow vector or array if necessary + + if (atom->nmax > nmax) { + nmax = atom->nmax; + if (nvalues == 1) { + memory->destroy(vector_atom); + memory->create(vector_atom,nmax,"rheo/property/atom:vector"); + } else { + memory->destroy(array_atom); + memory->create(array_atom,nmax,nvalues,"rheo/property/atom:array"); + } + } + + // fill vector or array with per-atom values + + if (nvalues == 1) { + buf = vector_atom; + (this->*pack_choice[0])(0); + } else { + if (nmax) buf = &array_atom[0][0]; + else buf = nullptr; + for (int n = 0; n < nvalues; n++) + (this->*pack_choice[n])(n); + } +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based array +------------------------------------------------------------------------- */ + +double ComputeRHEOPropertyAtom::memory_usage() +{ + double bytes = (double)nmax * nvalues * sizeof(double); + return bytes; +} + +/* ---------------------------------------------------------------------- + one method for every keyword compute rheo/property/atom can output + the atom property is packed into buf starting at n with stride nvalues + customize a new keyword by adding a method +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_phase(int n) +{ + int *status = atom->status; + 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); + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_chi(int n) +{ + double *chi = compute_interface->chi; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = chi[i]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_surface(int n) +{ + int *status = atom->status; + 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); + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_surface_r(int n) +{ + double *rsurface = compute_surface->rsurface; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = rsurface[i]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_surface_divr(int n) +{ + double *divr = compute_surface->divr; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = divr[i]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_surface_nx(int n) +{ + double **nsurface = compute_surface->nsurface; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = nsurface[i][0]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_surface_ny(int n) +{ + double **nsurface = compute_surface->nsurface; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = nsurface[i][1]; + else buf[n] = 0.0; + n += nvalues; + } +} + + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_surface_nz(int n) +{ + double **nsurface = compute_surface->nsurface; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = nsurface[i][2]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_coordination(int n) +{ + int *coordination = compute_kernel->coordination; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = coordination[i]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_cv(int n) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = fix_thermal->calc_cv(i); + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_shift_vx(int n) +{ + double **vshift = compute_vshift->vshift; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = vshift[i][0]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_shift_vy(int n) +{ + double **vshift = compute_vshift->vshift; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = vshift[i][1]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_shift_vz(int n) +{ + double **vshift = compute_vshift->vshift; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = vshift[i][2]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOPropertyAtom::pack_atom_style(int n) +{ + atom->avec->pack_property_atom(index[n],&buf[n],nvalues,groupbit); +} diff --git a/src/compute_rheo_property_atom.h b/src/compute_rheo_property_atom.h new file mode 100644 index 0000000000..accb7e9156 --- /dev/null +++ b/src/compute_rheo_property_atom.h @@ -0,0 +1,71 @@ +/* -*- 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/property/atom,ComputeRHEOPropertyAtom); +// clang-format on +#else + +#ifndef LMP_COMPUTE_RHEO_PROPERTY_ATOM_H +#define LMP_COMPUTE_RHEO_PROPERTY_ATOM_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeRHEOPropertyAtom : public Compute { + public: + ComputeRHEOPropertyAtom(class LAMMPS *, int, char **); + ~ComputeRHEOPropertyAtom() override; + void init() override; + void compute_peratom() override; + double memory_usage() override; + + private: + int nvalues, nmax; + int thermal_flag, interface_flag, surface_flag, shift_flag; + int *index; + double *buf; + + typedef void (ComputeRHEOPropertyAtom::*FnPtrPack)(int); + FnPtrPack *pack_choice; // ptrs to pack functions + + void pack_phase(int); + void pack_chi(int); + void pack_surface(int); + void pack_surface_r(int); + void pack_surface_divr(int); + void pack_surface_nx(int); + void pack_surface_ny(int); + void pack_surface_nz(int); + void pack_coordination(int); + void pack_cv(int); + void pack_shift_vx(int); + void pack_shift_vy(int); + void pack_shift_vz(int); + void pack_atom_style(int); + + class FixRHEO *fix_rheo; + class FixRHEOThermal *fix_thermal; + class ComputeRHEOInterface *compute_interface; + class ComputeRHEOKernel *compute_kernel; + class ComputeRHEOSurface *compute_surface; + class ComputeRHEOVShift *compute_vshift; + +}; + +} // namespace LAMMPS_NS + +#endif +#endif