From de0e4bb170f7a2f7ae59616ca4f4d203e6eccc68 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Thu, 20 Apr 2023 16:17:04 -0600 Subject: [PATCH] Rho sum compute --- src/RHEO/compute_rheo_kernel.cpp | 234 ++++++++++++++++-------------- src/RHEO/compute_rheo_kernel.h | 3 + src/RHEO/compute_rheo_rho_sum.cpp | 184 +++++++++++++++++++++++ src/RHEO/compute_rheo_rho_sum.h | 50 +++++++ src/RHEO/compute_rheo_surface.cpp | 2 +- src/RHEO/compute_rheo_surface.h | 17 +-- src/RHEO/fix_rheo.cpp | 6 +- 7 files changed, 371 insertions(+), 125 deletions(-) create mode 100644 src/RHEO/compute_rheo_rho_sum.cpp create mode 100644 src/RHEO/compute_rheo_rho_sum.h diff --git a/src/RHEO/compute_rheo_kernel.cpp b/src/RHEO/compute_rheo_kernel.cpp index 2fd189e1e4..c0b61daae3 100644 --- a/src/RHEO/compute_rheo_kernel.cpp +++ b/src/RHEO/compute_rheo_kernel.cpp @@ -47,21 +47,46 @@ using namespace MathExtra; enum {QUINTIC, CRK0, CRK1, CRK2}; #define DELTA 2000 -Todo: convert delx dely delz to an array -Should vshift be using kernel quintic? -Move away from h notation, use cut? - /* ---------------------------------------------------------------------- */ ComputeRHEOKernel::ComputeRHEOKernel(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), list(nullptr), C(nullptr), C0(nullptr), coordination(nullptr), compute_interface(nullptr) { - if (narg != 3) error->all(FLERR,"Illegal compute rheo/kernel command"); + if (narg != 4) error->all(FLERR,"Illegal compute rheo/kernel command"); + + kernel_style = (SubModelType) utils::inumeric(FLERR,arg[3],false,lmp); + + + if (kernel_style == FixRHEO::QUINTIC) { + correction_order = -1; + } else if (kernel_style == FixRHEO::CRK0) { + correction_order = 0; + } else if (kernel_style == FixRHEO::CRK1) { + correction_order = 1; + } else if (kernel_style == FixRHEO::CRK2) { + correction_order = 2; + } - comm_forward = 1; // Always minimum for coordination solid_flag = 0; dim = domain->dimension; + + comm_forward = 1; + ncor = 0; + Mdim = 0; + if (kernel_type == CRK1) { + Mdim = 1 + dim; + ncor = 1 + dim; + comm_forward = ncor * Mdim; + } else if (kernel_type == CRK2) { + //Polynomial basis size (up to quadratic order) + Mdim = 1 + dim + dim * (dim + 1) / 2; + //Number of sets of correction coefficients (1 x y xx yy) + z zz (3D) + ncor = 1 + 2 * dim; + comm_forward = ncor * Mdim; + } + + comm_forward_save = comm_forward; } /* ---------------------------------------------------------------------- */ @@ -93,22 +118,12 @@ void ComputeRHEOKernel::init() solid_flag = 1; } - kernel_style = fix_rheo->kernel_style; zmin = fix_rheo->zmin_kernel; h = fix_rheo->h; hsq = h * h; hinv = 1.0 / h; hsqinv = hinv * hinv; - if (kernel_style == FixRHEO::QUINTIC) { - correction_order = -1; - } else if (kernel_style == FixRHEO::CRK0) { - correction_order = 0; - } else if (kernel_style == FixRHEO::CRK1) { - correction_order = 1; - } else if (kernel_style == FixRHEO::CRK2) { - correction_order = 2; - } if (dim == 3) { pre_w = 0.002652582384864922 * 27.0 * ihsq * ih; @@ -133,23 +148,13 @@ void ComputeRHEOKernel::init() coordination = atom->ivector[index]; // Create local arrays for kernel arrays, I can't foresee a reason to print - comm_forward = 1; - ncor = 0; - Mdim = 0; + if (kernel_type == CRK0) { memory->create(C0, nmax_old, "rheo/kernel:C0"); } else if (kernel_type == CRK1) { - Mdim = 1 + dim; - ncor = 1 + dim; memory->create(C, nmax_old, ncor, Mdim, "rheo/kernel:C"); - comm_forward = ncor * Mdim; } else if (kernel_type == CRK2) { - //Polynomial basis size (up to quadratic order) - Mdim = 1 + dim + dim * (dim + 1) / 2; - //Number of sets of correction coefficients (1 x y xx yy) + z zz (3D) - ncor = 1 + 2 * dim; memory->create(C, nmax_old, ncor, Mdim, "rheo/kernel:C"); - comm_forward = ncor * Mdim; } } @@ -491,6 +496,8 @@ void ComputeRHEOKernel::compute_peratom() gsl_error_flag = 0; gsl_error_tags.clear(); + if (kernel_type == FixRHEO::QUINTIC) return; + int i, j, ii, jj, inum, jnum, g, a, b, gsl_error; double xtmp, ytmp, ztmp, r, rsq, w, vj; double dx[3]; @@ -513,48 +520,11 @@ void ComputeRHEOKernel::compute_peratom() firstneigh = list->firstneigh; // Grow arrays if necessary - int nmax = atom->nmax; - if (nmax_old < nmax) - memory->grow(coordination, nmax, "atom:rheo_coordination"); + if (nmax_old < atom->nmax) grow_arrays(atom->nmax); - if (kernel_type == FixRHEO::CRK0) { - memory->grow(C0, nmax, "rheo/kernel:C0"); - } else if (correction_order > 0) { - memory->grow(C, nmax, ncor, Mdim, "rheo/kernel:C"); - } - - nmax_old = atom->nmax; - } - - if (kernel_type == FixRHEO::QUINTIC) { - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - - jlist = firstneigh[i]; - jnum = numneigh[i]; - coordination[i] = 0; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - j &= NEIGHMASK; - - dx[0] = xtmp - x[j][0]; - dx[1] = ytmp - x[j][1]; - dx[2] = ztmp - x[j][2]; - rsq = lensq(dx); - - if (rsq < hsq) { - coordination[i] += 1; - } - } - } - } else if (kernel_type == FixRHEO::CRK0) { + if (kernel_type == FixRHEO::CRK0) { double M; - for (ii = 0; ii < inum; ii++) { i = ilist[ii]; xtmp = x[i][0]; @@ -566,7 +536,6 @@ void ComputeRHEOKernel::compute_peratom() //Initialize M to zero: M = 0; - coordination[i] = 0; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; @@ -584,7 +553,6 @@ void ComputeRHEOKernel::compute_peratom() vj = mass[type[j]] / compute_interface->correct_rho(j,i); } else vj = mass[type[j]] / rho[j]; - coordination[i] += 1; M += w * vj; } } @@ -613,7 +581,6 @@ void ComputeRHEOKernel::compute_peratom() M[a * Mdim + b] = 0; } } - coordination[i] = 0; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; @@ -658,8 +625,6 @@ void ComputeRHEOKernel::compute_peratom() } } - coordination[i] += 1; - // Populate the upper triangle for (a = 0; a < Mdim; a++) { for (b = a; b < Mdim; b++) { @@ -733,7 +698,74 @@ void ComputeRHEOKernel::compute_peratom() } // communicate calculated quantities - comm->forward_comm_compute(this); + comm_stage = 1; + comm_forward = comm_forward_save; + comm->forward_comm(this); +} + + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOKernel::compute_coordination() +{ + int i, j, ii, jj, inum, jnum; + double xtmp, ytmp, ztmp, rsq; + double dx[3]; + + double **x = atom->x; + + int *ilist, *jlist, *numneigh, **firstneigh; + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // Grow arrays if necessary + if (nmax_old < atom->nmax) grow_arrays(atom->nmax); + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + jlist = firstneigh[i]; + jnum = numneigh[i]; + coordination[i] = 0; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + dx[0] = xtmp - x[j][0]; + dx[1] = ytmp - x[j][1]; + dx[2] = ztmp - x[j][2]; + rsq = lensq(dx); + + if (rsq < hsq) + coordination[i] += 1; + } + } + + // communicate calculated quantities + comm_stage = 0; + comm_forward = 1; + comm->forward_comm(this); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEOKernel::grow_arrays(int nmax) +{ + memory->grow(coordination, nmax, "atom:rheo_coordination"); + + if (kernel_type == FixRHEO::CRK0) { + memory->grow(C0, nmax, "rheo/kernel:C0"); + } else if (correction_order > 0) { + memory->grow(C, nmax, ncor, Mdim, "rheo/kernel:C"); + } + + nmax_old = nmax; } /* ---------------------------------------------------------------------- */ @@ -743,26 +775,19 @@ int ComputeRHEOKernel::pack_forward_comm(int n, int *list, double *buf, { int i,j,k,m,a,b; m = 0; - if (correction_order > 0) { - for (i = 0; i < n; i++) { - j = list[i]; - for (a = 0; a < ncor; a++) { - for (b = 0; b < Mdim; b++) { - buf[m++] = C[j][a][b]; - } + + for (i = 0; i < n; i++) { + j = list[i]; + if (comm_stage == 0) { + buf[m++] = coordination[j]; + } else { + if (kernel_type == FixRHEO::CRK0) { + buf[m++] = C0[j]; + } else { + for (a = 0; a < ncor; a++) + for (b = 0; b < Mdim; b++) + buf[m++] = C[j][a][b]; } - buf[m++] = coordination[j]; - } - } else if (kernel_type == FixRHEO::CRK0) { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = C0[j]; - buf[m++] = coordination[j]; - } - } else { - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = coordination[j]; } } return m; @@ -775,23 +800,18 @@ void ComputeRHEOKernel::unpack_forward_comm(int n, int first, double *buf) int i, k, m, last,a,b; m = 0; last = first + n; - if (correction_order > 0) { - for (i = first; i < last; i++) { - for (a = 0; a < ncor; a++) { - for (b = 0; b < Mdim; b++) { - C[i][a][b] = buf[m++]; - } + + for (i = first; i < last; i++) { + if (comm_stage == 0) { + coordination[i] = buf[m++]; + } else { + if (kernel_type == FixRHEO::CRK0) { + C0[i] = buf[m++]; + } else { + for (a = 0; a < ncor; a++) + for (b = 0; b < Mdim; b++) + C[i][a][b] = buf[m++]; } - coordination[i] = buf[m++]; - } - } else if (kernel_type == FixRHEO::CRK0) { - for (i = first; i < last; i++) { - C0[i] = buf[m++]; - coordination[i] = buf[m++]; - } - } else { - for (i = first; i < last; i++) { - coordination[i] = buf[m++]; } } } diff --git a/src/RHEO/compute_rheo_kernel.h b/src/RHEO/compute_rheo_kernel.h index ed2b6ff5f2..4fbdb966b4 100644 --- a/src/RHEO/compute_rheo_kernel.h +++ b/src/RHEO/compute_rheo_kernel.h @@ -35,16 +35,19 @@ class ComputeRHEOKernel : public Compute { int pack_forward_comm(int, int *, double *, int, int *) override; void unpack_forward_comm(int, int, double *) override; double memory_usage() override; + void compute_coordination(); double calc_w(int,int,double,double,double,double); double calc_dw(int,int,double,double,double,double); double calc_w_quintic(int,int,double,double,double,double); double calc_dw_quintic(int,int,double,double,double,double,double *,double *); + void grow_arrays(int); double dWij[3], dWji[3], Wij, Wji; int correction_order; int *coordination; private: + int comm_stage, comm_forward_save; int solid_flag; int gsl_error_flag; std::unordered_set gsl_error_tags; diff --git a/src/RHEO/compute_rheo_rho_sum.cpp b/src/RHEO/compute_rheo_rho_sum.cpp new file mode 100644 index 0000000000..fafd948538 --- /dev/null +++ b/src/RHEO/compute_rheo_rho_sum.cpp @@ -0,0 +1,184 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + 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_rho_sum.h" + +#include "atom.h" +#include "comm.h" +#include "compute_rheo_kernel.h" +#include "error.h" +#include "fix_rheo.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputeRHEORhoSum::ComputeRHEORhoSum(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), fix_rheo(nullptr), compute_kernel(nullptr) +{ + if (narg != 3) error->all(FLERR,"Illegal compute RHEO/rho command"); + + comm_forward = 1; + comm_reverse = 1; +} + +/* ---------------------------------------------------------------------- */ + +ComputeRHEORhoSum::~ComputeRHEORhoSum() {} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEORhoSum::init() +{ + compute_kernel = fix_rheo->compute_kernel; + cut = fix_rheo->cut; + cutsq = cut * cut; + + // need an occasional half neighbor list + neighbor->add_request(this, NeighConst::REQ_HALF); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEORhoSum::init_list(int /*id*/, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- */ + + +void ComputeRHEORhoSum::compute_peratom() +{ + int i, j, ii, jj, inum, jnum, itype, jtype; + double xtmp, ytmp, ztmp, delx, dely, delz; + int *ilist, *jlist, *numneigh, **firstneigh; + double rsq, w; + + int nlocal = atom->nlocal; + + double **x = atom->x; + double *rho = atom->rho; + int *type = atom->type; + int *status = atom->status; + double *mass = atom->mass; + int newton = force->newton; + + double jmass; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + int nall = nlocal + atom->nghost; + + // initialize arrays, local with quintic self-contribution, ghosts are zeroed + for (i = 0; i < nlocal; i++) { + w = compute_kernel->calc_w_quintic(i, i, 0.0, 0.0, 0.0, 0.0); + rho[i] += w * mass[type[i]]; + } + + for (i = nlocal; i < nall; i++) rho[i] = 0.0; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx * delx + dely * dely + delz * delz; + if (rsq < cutsq) { + w = compute_kernel->calc_w(i, j, delx, dely, delz, sqrt(rsq)); + rho[i] += w * mass[type[i]]; + if (newton || j < nlocal) { + rho[j] += w * mass[type[j]]; + } + } + } + } + + if (newton) comm->reverse_comm(this); + comm->forward_comm(this); +} + +/* ---------------------------------------------------------------------- */ + +int ComputeRHEORhoSum::pack_forward_comm(int n, int *list, double *buf, + int /*pbc_flag*/, int * /*pbc*/) +{ + int i,j,k,m; + double * rho = atom->rho; + m = 0; + + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = coordination[j]; + } + return m; +} + +/* ---------------------------------------------------------------------- */ +void ComputeRHEORhoSum::unpack_forward_comm(int n, int first, double *buf) +{ + int i, k, m, last; + double * rho = atom->rho; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + rho[i] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + +int ComputeRHEORhoSum::pack_reverse_comm(int n, int first, double *buf) +{ + int i,k,m,last; + double *rho = atom->rho; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + buf[m++] = rho[i]; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRHEORhoSum::unpack_reverse_comm(int n, int *list, double *buf) +{ + int i,k,j,m; + double *rho = atom->rho; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + rho[j] += buf[m++]; + } +} diff --git a/src/RHEO/compute_rheo_rho_sum.h b/src/RHEO/compute_rheo_rho_sum.h new file mode 100644 index 0000000000..a411d5ed29 --- /dev/null +++ b/src/RHEO/compute_rheo_rho_sum.h @@ -0,0 +1,50 @@ +/* -*- 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/RHO/SUM,ComputeRHEORhoSum) +// clang-format on +#else + +#ifndef LMP_COMPUTE_RHEO_RHO_SUM_H +#define LMP_COMPUTE_RHEO_RHO_SUM_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeRHEORhoSum : public Compute { + public: + ComputeRHEORhoSum(class LAMMPS *, int, char **); + ~ComputeRHEORhoSum() override; + void init() override; + void init_list(int, class NeighList *) override; + void compute_peratom() override; + int pack_forward_comm(int, int *, double *, int, int *) override; + void unpack_forward_comm(int, int, double *) override; + int pack_reverse_comm(int, int, double *) override; + void unpack_reverse_comm(int, int *, double *) override; + + private: + double cut, cutsq; + + class NeighList *list; + class FixRHEO *fix_rheo; + class ComputeRHEOKernel *compute_kernel; +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/RHEO/compute_rheo_surface.cpp b/src/RHEO/compute_rheo_surface.cpp index 5b3aeabf26..9a94ea6a76 100644 --- a/src/RHEO/compute_rheo_surface.cpp +++ b/src/RHEO/compute_rheo_surface.cpp @@ -36,7 +36,7 @@ using namespace LAMMPS_NS; using namespace FixConst; using namespace MathExtra; -#define epsilon 1e-10; +#define EPSILON 1e-10; /* ---------------------------------------------------------------------- */ diff --git a/src/RHEO/compute_rheo_surface.h b/src/RHEO/compute_rheo_surface.h index 8d0b681c6f..00cfb56b31 100644 --- a/src/RHEO/compute_rheo_surface.h +++ b/src/RHEO/compute_rheo_surface.h @@ -27,7 +27,7 @@ namespace LAMMPS_NS { class ComputeRHEOSurface : public Compute { public: ComputeRHEOSurface(class LAMMPS *, int, char **); - ~ComputeRHEOSurface(); + ~ComputeRHEOSurface() override; void init() override; void init_list(int, class NeighList *) override; void compute_peratom() override; @@ -44,26 +44,13 @@ class ComputeRHEOSurface : public Compute { double **B, **gradC, *divr; int threshold_style, comm_stage; - double divR_limit; - int coord_limit; - class NeighList *list; class FixRHEO *fix_rheo; class ComputeRHEOKernel *compute_kernel; class ComputeRHEOSolids *compute_solids; }; -} +} // namespace LAMMPS_NS #endif #endif - -/* ERROR/WARNING messages: - -E: Illegal ... command - -Self-explanatory. Check the input script syntax and compare to the -documentation for the command. You can use -echo screen as a -command-line option when running LAMMPS to see the offending line. - -*/ diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index 5358cc4cba..e2f9467b34 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -141,7 +141,7 @@ FixRHEO::~FixRHEO() void FixRHEO::post_constructor() { - compute_kernel = dynamic_cast(modify->add_compute("rheo_kernel all RHEO/KERNEL")); + compute_kernel = dynamic_cast(modify->add_compute("rheo_kernel all RHEO/KERNEL {}", kernel_style)); compute_kernel->fix_rheo = this; std::string cmd = "rheo_grad all RHEO/GRAD velocity rho viscosity"; @@ -150,7 +150,7 @@ void FixRHEO::post_constructor() compute_grad->fix_rheo = this; if (rhosum_flag) { - compute_rhosum = dynamic_cast(modify->add_compute("rheo_rhosum all RHEO/RHO/SUM")); + compute_rhosum = dynamic_cast(modify->add_compute("rheo_rho_sum all RHEO/RHO/SUM")); compute_rhosum->fix_rheo = this; } @@ -364,6 +364,8 @@ void FixRHEO::initial_integrate(int /*vflag*/) void FixRHEO::pre_force(int /*vflag*/) { + compute_kernel->compute_coordination(); // Needed for rho sum + if (rhosum_flag) compute_rhosum->compute_peratom();