diff --git a/src/RHEO/compute_rheo_grad.cpp b/src/RHEO/compute_rheo_grad.cpp new file mode 100644 index 0000000000..33dd3dc3bc --- /dev/null +++ b/src/RHEO/compute_rheo_grad.cpp @@ -0,0 +1,457 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "compute_rheo_grad.h" + +#include "atom.h" +#include "comm.h" +#include "compute_rheo_kernel.h" +#include "compute_rheo_solids.h" +#include "domain.h" +#include "error.h" +#include "fix_rheo.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "pair.h" +#include "memory.h" +#include "modify.h" +#include "update.h" +#include "utils.h" + +#include +#include + +using namespace LAMMPS_NS; +enum{COMMGRAD, COMMFIELD}; + +/* ---------------------------------------------------------------------- */ + +ComputeRHEOGrad::ComputeRHEOGrad(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), compute_interface(nullptr), compute_kernel(nullptr) +{ + if (narg < 4) error->all(FLERR,"Illegal compute rheo/grad command"); + + + velocity_flag = temperature_flag = rho_flag = eta_flag = 0; + for (int iarg = 3; iarg < narg; iarg ++) { + if (strcmp(arg[iarg],"velocity") == 0) velocity_flag = 1; + else if (strcmp(arg[iarg],"rho") == 0) rho_flag = 1; + else if (strcmp(arg[iarg],"temperature") == 0) temperature_flag = 1; + else if (strcmp(arg[iarg],"viscosity") == 0) eta_flag = 1; + else error->all(FLERR, "Illegal compute rheo/grad command, {}", arg[iarg]); + } + + dim = domain->dimension; + + ncomm_grad = 0; + ncomm_field = 0; + comm_reverse = 0; + + std::string fix_cmd = "rheo_grad_property_atom all property/atom" + + if (velocity_flag) { + ncomm_grad += dim * dim; + ncomm_field += dim; + comm_reverse += dim * dim; + fix_cmd += " d2_gradv 9" + } + + if (rho_flag) { + ncomm_grad += dim; + ncomm_field += 1; + comm_reverse += dim; + fix_cmd += " d2_gradr 3" + } + + if (temperature_flag) { + ncomm_grad += dim; + ncomm_field += 1; + comm_reverse += dim; + fix_cmd += " d2_gradt 3" + } + + if (eta_flag) { + ncomm_grad += dim; + comm_reverse += dim; + fix_cmd += " d2_gradn 3" + } + + comm_forward = ncomm_grad; + + modify->add_fix(fix_cmd); + + int tmp1, tmp2, index; + if (velocity_flag) { + index = atom->find_custom("gradv", tmp1, tmp2); + gradv = atom->darray[index]; + } + + if (rho_flag) { + index = atom->find_custom("gradr", tmp1, tmp2); + gradr = atom->darray[index]; + } + + if (temperature_flag) { + index = atom->find_custom("gradt", tmp1, tmp2); + gradt = atom->darray[index]; + } + + if (eta_flag) { + index = atom->find_custom("gradn", tmp1, tmp2); + gradn = atom->darray[index]; + } +} + +/* ---------------------------------------------------------------------- */ + +ComputeRHEOGrad::~ComputeRHEOGrad() +{ + modify->delete_fix("rheo_grad_property_atom"); +} + + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOGrad::init() +{ + neighbor->add_request(this, NeighConst::REQ_DEFAULT); + + cut = fix_rheo->cut; + cutsq = cut * cut; + rho0 = fix_rheo->rho0; + compute_kernel = fix_rheo->compute_kernel; + compute_interface = fix_rheo->compute_interface; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOGrad::init_list(int /*id*/, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOGrad::compute_peratom() +{ + int i, j, k, ii, jj, jnum, itype, jtype, a, b; + double xtmp, ytmp, ztmp, delx, dely, delz; + double rsq, imass, jmass; + double rhoi, rhoj, Voli, Volj, drho, dT, deta; + double vij[3]; + double wp, *dWij, *dWji; + + int inum, *ilist, *numneigh, **firstneigh; + int *jlist; + int nlocal = atom->nlocal; + + double **x = atom->x; + double **v = atom->v; + double *rho = atom->rho; + double *temperature = atom->temperature; + double *eta = atom->viscosity; + int *status = atom->status; + int *type = atom->type; + double *mass = atom->mass; + int newton = force->newton; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // initialize arrays + for (i = 0; i < nmax; i++) { + if (velocity_flag) { + for (k = 0; k < dim * dim; k++) + gradv[i][k] = 0.0; + } + if (rho_flag) { + for (k = 0; k < dim; k++) + gradr[i][k] = 0.0; + } + if (temperature_flag) { + for (k = 0; k < dim; k++) + gradt[i][k] = 0.0; + } + if (eta_flag) { + for (k = 0; k < dim; k++) + gradn[i][k] = 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]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + jtype = type[j]; + 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) { + rhoi = rho[i]; + rhoj = rho[j]; + + // Add corrections for walls + if ((status[i] & FixRHEO::STATUS_FLUID) && !(status[j] & FixRHEO::STATUS_FLUID)) { + compute_interface->correct_v(v[i], v[j], vi, i, j); + rhoj = compute_interface->correct_rho(j,i); + } else if (!(status[i] & FixRHEO::STATUS_FLUID) && (status[j] & FixRHEO::STATUS_FLUID)) { + compute_interface->correct_v(v[j], v[i], vj, j, i); + rhoi = compute_interface->correct_rho(i,j); + } else if (!(status[i] & FixRHEO::STATUS_FLUID) && !(status[j] & FixRHEO::STATUS_FLUID)) { + rhoi = rho0; + rhoj = rho0; + } + + Voli = mass[itype] / rhoi; + Volj = mass[jtype] / rhoj; + + vij[0] = v[i][0] - v[j][0]; + vij[1] = v[i][1] - v[j][1]; + vij[2] = v[i][2] - v[j][2]; + + if (rho_flag) drho = rhoi - rhoj; + if (temperature_flag) dT = temperature[i] - temperature[j]; + if (eta_flag) deta = eta[i] - eta[j]; + + wp = compute_kernel->calc_dw(i, j, delx, dely, delz, sqrt(rsq)); + dWij = compute_kernel->dWij; + dWji = compute_kernel->dWji; + + for (a = 0; a < dim; a++) { + for (b = 0; b < dim; b++) { + if (velocity_flag) // uxx uxy uxz uyx uyy uyz uzx uzy uzz + gradv[i][a * dim + b] -= vij[a] * Volj * dWij[b]; + } + + if (rho_flag) // P,x P,y P,z + gradr[i][a] -= drho * Volj * dWij[a]; + + if (temperature_flag) // T,x T,y T,z + gradt[i][a] -= dT * Volj * dWij[a]; + + if (eta_flag) // n,x n,y n,z + gradn[i][a] -= deta * Volj * dWij[a]; + } + + if (newton || j < nlocal) { + for (a = 0; a < dim; a++) { + for (b = 0; b < dim; b++) { + if (velocity_flag) // uxx uxy uxz uyx uyy uyz uzx uzy uzz + gradv[j][a * dim + b] += vij[a] * Voli * dWji[b]; + } + + if (rho_flag) // P,x P,y P,z + gradr[j][a] += drho * Voli * dWji[a]; + + if (temperature_flag) // T,x T,y T,z + gradt[j][a] += dT * Voli * dWji[a]; + + if (eta_flag) // n,x n,y n,z + gradn[j][a] += deta * Voli * dWji[a]; + } + } + } + } + } + + if (newton) comm->reverse_comm(this); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOGrad::forward_gradients() +{ + comm_stage = COMMGRAD; + comm_forward = ncomm_grad; + comm->forward_comm(this); +} + + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOGrad::forward_fields() +{ + comm_stage = COMMFIELD; + comm_forward = ncomm_field; + comm->forward_comm(this); +} + +/* ---------------------------------------------------------------------- */ + +int ComputeRHEOGrad::pack_forward_comm(int n, int *list, double *buf, + int /*pbc_flag*/, int * /*pbc*/) +{ + int i,j,k,m; + double *rho = atom->rho; + double *temperature = atom->temperature; + double *eta = atom->viscosity; + double **v = atom->v; + + m = 0; + + for (i = 0; i < n; i++) { + j = list[i]; + if (comm_stage == COMMGRAD) { + + if (velocity_flag){ + for (k = 0; k < dim * dim; k++) + buf[m++] = gradv[j][k]; + } + + if (rho_flag) { + for (k = 0; k < dim; k++) + buf[m++] = gradr[j][k]; + } + + if (temperature_flag) { + for (k = 0; k < dim; k++) + buf[m++] = gradt[j][k]; + } + + if (eta_flag){ + for (k = 0; k < dim; k++) + buf[m++] = gradn[j][k]; + } + } else if (comm_stage == COMMFIELD) { + + if (velocity_flag) { + for (k = 0; k < dim; k++) + buf[m++] = v[j][k]; + } + + if (rho_flag) { + buf[m++] = rho[j]; + } + + if (temperature_flag) { + buf[m++] = temperature[j]; + } + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOGrad::unpack_forward_comm(int n, int first, double *buf) +{ + int i, k, m, last; + double * rho = atom->rho; + double * temperature = atom->temperature; + double ** v = atom->v; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + if (comm_stage == COMMGRAD) { + if (velocity_flag) { + for (k = 0; k < dim * dim; k++) + gradv[i][k] = buf[m++]; + } + if (rho_flag) { + for (k = 0; k < dim; k++) + gradr[i][k] = buf[m++]; + } + if (temperature_flag) { + for (k = 0; k < dim; k++) + gradt[i][k] = buf[m++]; + } + if (eta_flag) { + for (k = 0; k < dim; k++) + gradn[i][k] = buf[m++]; + } + } else if (comm_stage == COMMFIELD) { + if (velocity_flag) { + for (k = 0; k < dim; k++) + v[i][k] = buf[m++]; + } + if (rho_flag) { + rho[i] = buf[m++]; + } + if (temperature_flag) { + temperature[i] = buf[m++]; + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +int ComputeRHEOGrad::pack_reverse_comm(int n, int first, double *buf) +{ + int i,k,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + if (velocity_flag) { + for (k = 0; k < dim * dim; k++) + buf[m++] = gradv[i][k]; + } + if (rho_flag) { + for (k = 0; k < dim; k++) + buf[m++] = gradr[i][k]; + } + if (temperature_flag) { + for (k = 0; k < dim; k++) + buf[m++] = gradt[i][k]; + } + if (eta_flag) { + for (k = 0; k < dim; k++) + buf[m++] = gradn[i][k]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOGrad::unpack_reverse_comm(int n, int *list, double *buf) +{ + int i,k,j,m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + if (velocity_flag) { + for (k = 0; k < dim * dim; k++) + gradv[j][k] += buf[m++]; + } + if (rho_flag) { + for (k = 0; k < dim; k++) + gradr[j][k] += buf[m++]; + } + if (temperature_flag) { + for (k = 0; k < dim; k++) + gradt[j][k] += buf[m++]; + } + if (eta_flag) { + for (k = 0; k < dim; k++) + gradn[j][k] += buf[m++]; + } + } +} diff --git a/src/RHEO/compute_rheo_grad.h b/src/RHEO/compute_rheo_grad.h new file mode 100644 index 0000000000..8b8e28fd98 --- /dev/null +++ b/src/RHEO/compute_rheo_grad.h @@ -0,0 +1,62 @@ +/* -*- 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/grad,ComputeRHEOGrad) +// clang-format on +#else + +#ifndef LMP_COMPUTE_RHEO_GRAD_H +#define LMP_COMPUTE_RHEO_GRAD_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeRHEOGrad : public Compute { + public: + ComputeRHEOGrad(class LAMMPS *, int, char **); + ~ComputeRHEOGrad(); + 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 forward_gradients(); + void forward_fields(); + double **gradv; + double **gradr; + double **gradt; + double **gradn; + int stage; + + private: + int dim, comm_stage; + int ncomm_grad, ncomm_field; + double cut, cutsq, rho0; + class NeighList *list; + + class FixRHEO *fix_rheo; + class ComputeRHEOKernel *compute_kernel; + class ComputeRHEOInterface *compute_interface; + + int velocity_flag, temperature_flag, rho_flag, eta_flag; +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index fd476d4c9e..17108b4c91 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -114,6 +114,7 @@ void FixRHEO::post_constructor() compute_grad = dynamic_cast(modify->add_compute(fmt::format("rheo_grad all rheo/grad {} velocity rho viscosity temprature", cut))); else compute_grad = dynamic_cast(modify->add_compute(fmt::format("rheo_grad all rheo/grad {} velocity rho viscosity", cut))); + compute_grad->fix_rheo = this; compute_interface = dynamic_cast(modify->add_compute(fmt::format("rheo_interface all rheo/interface {}", cut))); diff --git a/src/fix_bond_history.cpp b/src/fix_bond_history.cpp index cae9dc744d..277df75085 100644 --- a/src/fix_bond_history.cpp +++ b/src/fix_bond_history.cpp @@ -97,7 +97,7 @@ void FixBondHistory::post_constructor() void FixBondHistory::update_atom_value(int i, int m, int idata, double value) { - if (idata >= ndata || m > nbond) error->all(FLERR, "Index exceeded in fix bond history"); + if (idata >= ndata || m > nbond) error->one(FLERR, "Index exceeded in fix bond history"); atom->darray[index][i][m * ndata + idata] = value; } @@ -105,7 +105,7 @@ void FixBondHistory::update_atom_value(int i, int m, int idata, double value) double FixBondHistory::get_atom_value(int i, int m, int idata) { - if (idata >= ndata || m > nbond) error->all(FLERR, "Index exceeded in fix bond history"); + if (idata >= ndata || m > nbond) error->one(FLERR, "Index exceeded in fix bond history"); return atom->darray[index][i][m * ndata + idata]; }