ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Rho sum compute

Browse Source
This commit is contained in:
jtclemm
2023-04-20 16:17:04 -06:00
parent a4d971df52
commit de0e4bb170
7 changed files with 371 additions and 125 deletions
Download Patch File Download Diff File Expand all files Collapse all files

222
src/RHEO/compute_rheo_kernel.cpp
View File

@ -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) {
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,27 +775,20 @@ 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++) {
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,24 +800,19 @@ 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++) {
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++];
}
}
}

3
src/RHEO/compute_rheo_kernel.h
View File

@ -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<tagint> gsl_error_tags;

184
src/RHEO/compute_rheo_rho_sum.cpp Normal file
View File

@ -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++];
}
}

50
src/RHEO/compute_rheo_rho_sum.h Normal file
View File

@ -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

2
src/RHEO/compute_rheo_surface.cpp
View File

@ -36,7 +36,7 @@ using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathExtra;
#define epsilon 1e-10;
#define EPSILON 1e-10;
/* ---------------------------------------------------------------------- */

17
src/RHEO/compute_rheo_surface.h
View File

@ -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.
*/

6
src/RHEO/fix_rheo.cpp
View File

@ -141,7 +141,7 @@ FixRHEO::~FixRHEO()
void FixRHEO::post_constructor()
{
compute_kernel = dynamic_cast<ComputeRHEOKernel *>(modify->add_compute("rheo_kernel all RHEO/KERNEL"));
compute_kernel = dynamic_cast<ComputeRHEOKernel *>(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<ComputeRHEORhoSum *>(modify->add_compute("rheo_rhosum all RHEO/RHO/SUM"));
compute_rhosum = dynamic_cast<ComputeRHEORhoSum *>(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();