Fix compiler/memory errors in tension, update properties in surface
This commit is contained in:
jtclemm
4
src/.gitignore
vendored
4
src/.gitignore
vendored
@ -219,6 +219,8 @@
|
||||
/fix_rheo_pressure.h
|
||||
/fix_rheo_stress.cpp
|
||||
/fix_rheo_stress.h
|
||||
/fix_rheo_tension.cpp
|
||||
/fix_rheo_tension.h
|
||||
/fix_rheo_thermal.cpp
|
||||
/fix_rheo_thermal.h
|
||||
/fix_rheo_viscosity.cpp
|
||||
@ -227,8 +229,6 @@
|
||||
/pair_rheo.h
|
||||
/pair_rheo_react.cpp
|
||||
/pair_rheo_react.h
|
||||
/pair_rheo_tension.cpp
|
||||
/pair_rheo_tension.h
|
||||
|
||||
/compute_grid.cpp
|
||||
/compute_grid.h
|
||||
|
||||
@ -355,8 +355,6 @@ int ComputeRHEOGrad::pack_forward_comm(int n, int *list, double *buf,
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
if (rho_flag)
|
||||
buf[m++] = rho[j];
|
||||
|
||||
|
||||
@ -87,27 +87,26 @@ void ComputeRHEOSurface::init()
|
||||
|
||||
cutsq = cut * cut;
|
||||
|
||||
// Create rsurface, divr, nsurface arrays if they don'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
|
||||
// Create rsurface, divr, nsurface arrays as custom atom properties,
|
||||
// can print with compute property/atom
|
||||
// no grow callback as there's no reason to copy/exchange data, manually grow
|
||||
// For B and gradC, create a local array since they are unlikely to be printed
|
||||
|
||||
int dim = domain->dimension;
|
||||
int tmp1, tmp2;
|
||||
int index = atom->find_custom("rheo_divr", tmp1, tmp2);
|
||||
if (index == -1) index = atom->add_custom("rheo_divr", 1, 0);
|
||||
divr = atom->dvector[index];
|
||||
index_divr = atom->find_custom("rheo_divr", tmp1, tmp2);
|
||||
if (index_divr == -1) index_divr = atom->add_custom("rheo_divr", 1, 0);
|
||||
divr = atom->dvector[index_divr];
|
||||
|
||||
index = atom->find_custom("rheo_rsurface", tmp1, tmp2);
|
||||
if (index == -1) index = atom->add_custom("rheo_rsurface", 1, 0);
|
||||
rsurface = atom->dvector[index];
|
||||
index_rsurf = atom->find_custom("rheo_rsurface", tmp1, tmp2);
|
||||
if (index_rsurf == -1) index_rsurf = atom->add_custom("rheo_rsurface", 1, 0);
|
||||
rsurface = atom->dvector[index_rsurf];
|
||||
|
||||
index = atom->find_custom("rheo_nsurface", tmp1, tmp2);
|
||||
if (index == -1) index = atom->add_custom("rheo_nsurface", 1, 3);
|
||||
nsurface = atom->darray[index];
|
||||
index_nsurf = atom->find_custom("rheo_nsurface", tmp1, tmp2);
|
||||
if (index_nsurf == -1) index_nsurf = atom->add_custom("rheo_nsurface", 1, dim);
|
||||
nsurface = atom->darray[index_nsurf];
|
||||
|
||||
nmax_store = atom->nmax;
|
||||
int dim = domain->dimension;
|
||||
memory->create(B, nmax_store, dim * dim, "rheo/surface:B");
|
||||
memory->create(gradC, nmax_store, dim * dim, "rheo/surface:gradC");
|
||||
|
||||
@ -148,32 +147,21 @@ void ComputeRHEOSurface::compute_peratom()
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
int nmax = atom->nmax;
|
||||
if (nmax_store <= nmax) {
|
||||
memory->grow(divr, nmax, "atom:rheo_divr");
|
||||
memory->grow(rsurface, nmax, "atom:rheo_rsurface");
|
||||
memory->grow(nsurface, nmax, 3, "atom:rheo_nsurface");
|
||||
// Grow and zero arrays
|
||||
if (nmax_store <= atom->nmax)
|
||||
grow_arrays(atom->nmax);
|
||||
|
||||
memory->grow(B, nmax, dim * dim, "rheo/surface:B");
|
||||
memory->grow(gradC, nmax, dim * dim, "rheo/surface:gradC");
|
||||
|
||||
nmax_store = atom->nmax;
|
||||
}
|
||||
size_t nbytes = nmax_store * sizeof(double);
|
||||
memset(&divr, 0, nbytes);
|
||||
memset(&rsurface, 0, nbytes);
|
||||
memset(&nsurface, 0, 3 * nbytes);
|
||||
memset(&gradC, 0, 3 * 3 * nbytes);
|
||||
memset(&B, 0, 3 * 3 * nbytes);
|
||||
|
||||
// Remove surface settings
|
||||
int nall = nlocal + atom->nghost;
|
||||
for (i = 0; i < nall; i++) {
|
||||
for (a = 0; a < dim; a++) {
|
||||
for (b = 0; b < dim; b++) {
|
||||
B[i][a * dim + b] = 0.0;
|
||||
gradC[i][a * dim + b] = 0.0;
|
||||
}
|
||||
nsurface[i][a] = 0.0;
|
||||
}
|
||||
divr[i] = 0.0;
|
||||
|
||||
// Remove surface settings
|
||||
for (i = 0; i < nall; i++)
|
||||
status[i] &= SURFACEMASK;
|
||||
}
|
||||
|
||||
// loop over neighbors to calculate the average orientation of neighbors
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
@ -424,3 +412,25 @@ void ComputeRHEOSurface::unpack_forward_comm(int n, int first, double *buf)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeRHEOSurface::grow_arrays(int nmax)
|
||||
{
|
||||
int dim = domain->dimension;
|
||||
|
||||
// Grow atom variables and reassign pointers
|
||||
memory->grow(atom->dvector[index_divr], nmax, "atom:rheo_divr");
|
||||
memory->grow(atom->dvector[index_rsurf], nmax, "atom:rheo_rsurface");
|
||||
memory->grow(atom->darray[index_nsurf], nmax, dim, "atom:rheo_nsurface");
|
||||
|
||||
divr = atom->dvector[index_divr];
|
||||
rsurface = atom->dvector[index_rsurf];
|
||||
nsurface = atom->darray[index_nsurf];
|
||||
|
||||
// Grow local variables
|
||||
memory->grow(B, nmax, dim * dim, "rheo/surface:B");
|
||||
memory->grow(gradC, nmax, dim * dim, "rheo/surface:gradC");
|
||||
|
||||
nmax_store = atom->nmax;
|
||||
}
|
||||
|
||||
@ -40,14 +40,18 @@ class ComputeRHEOSurface : public Compute {
|
||||
class FixRHEO *fix_rheo;
|
||||
|
||||
private:
|
||||
double cut, cutsq, rho0, threshold_divr;
|
||||
int surface_style, nmax_store, threshold_z, threshold_splash, interface_flag;
|
||||
double **B, **gradC;
|
||||
int threshold_style, comm_stage;
|
||||
int index_divr, index_rsurf, index_nsurf;
|
||||
|
||||
double cut, cutsq, rho0, threshold_divr;
|
||||
double **B, **gradC;
|
||||
|
||||
class NeighList *list;
|
||||
class ComputeRHEOKernel *compute_kernel;
|
||||
class ComputeRHEOInterface *compute_interface;
|
||||
|
||||
void grow_arrays(int);
|
||||
};
|
||||
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
452
src/RHEO/fix_rheo_tension.cpp
Normal file
452
src/RHEO/fix_rheo_tension.cpp
Normal file
@ -0,0 +1,452 @@
|
||||
/* ----------------------------------------------------------------------
|
||||
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 "fix_rheo_tension.h"
|
||||
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "compute_rheo_kernel.h"
|
||||
#include "compute_rheo_interface.h"
|
||||
#include "domain.h"
|
||||
#include "error.h"
|
||||
#include "fix_rheo.h"
|
||||
#include "force.h"
|
||||
#include "math_extra.h"
|
||||
#include "memory.h"
|
||||
#include "modify.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "update.h"
|
||||
#include "utils.h"
|
||||
|
||||
#include <cmath>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace RHEO_NS;
|
||||
using namespace MathExtra;
|
||||
using namespace FixConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixRHEOTension::FixRHEOTension(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg), compute_kernel(nullptr), compute_interface(nullptr), fix_rheo(nullptr)
|
||||
{
|
||||
if (narg != 4) error->all(FLERR,"Illegal fix command");
|
||||
alpha = utils::numeric(FLERR,arg[3],false,lmp);
|
||||
|
||||
comm_forward = 3;
|
||||
comm_reverse = 3;
|
||||
|
||||
// Create cgrad, n, and divr arrays as custom atom properties,
|
||||
// can print with compute property/atom
|
||||
// no grow callback as there's no reason to copy/exchange data, manually grow
|
||||
// For norm, create a local array since they are unlikely to be printed
|
||||
|
||||
int tmp1, tmp2;
|
||||
index_cgradt = atom->find_custom("cgrad_rheo_tension", tmp1, tmp2);
|
||||
if (index_cgradt == -1) index_cgradt = atom->add_custom("cgrad_rheo_tension", 1, 3);
|
||||
cgradt = atom->darray[index_cgradt];
|
||||
|
||||
index_nt = atom->find_custom("n_rheo_tension", tmp1, tmp2);
|
||||
if (index_nt == -1) index_nt = atom->add_custom("n_rheo_tension", 1, 3);
|
||||
nt = atom->darray[index_nt];
|
||||
|
||||
index_divnt = atom->find_custom("divn_rheo_tension", tmp1, tmp2);
|
||||
if (index_divnt == -1) index_divnt = atom->add_custom("divn_rheo_tension", 1, 0);
|
||||
divnt = atom->dvector[index_divnt];
|
||||
|
||||
norm = nullptr;
|
||||
nmax_store = 0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixRHEOTension::~FixRHEOTension()
|
||||
{
|
||||
// Remove custom property if it exists
|
||||
int tmp1, tmp2, index;
|
||||
|
||||
index = atom->find_custom("cgrad_rheo_tension", tmp1, tmp2);
|
||||
if (index != -1) atom->remove_custom(index, 1, 3);
|
||||
|
||||
index = atom->find_custom("n_rheo_tension", tmp1, tmp2);
|
||||
if (index != -1) atom->remove_custom(index, 1, 3);
|
||||
|
||||
index = atom->find_custom("divn_rheo_tension", tmp1, tmp2);
|
||||
if (index != -1) atom->remove_custom(index, 1, 0);
|
||||
|
||||
memory->destroy(norm);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixRHEOTension::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= POST_FORCE;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixRHEOTension::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/tension");
|
||||
fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
|
||||
|
||||
compute_kernel = fix_rheo->compute_kernel;
|
||||
compute_interface = fix_rheo->compute_interface;
|
||||
interface_flag = fix_rheo->interface_flag;
|
||||
h = fix_rheo->h;
|
||||
rho0 = fix_rheo->rho0;
|
||||
|
||||
hsq = h * h;
|
||||
|
||||
neighbor->add_request(this, NeighConst::REQ_DEFAULT);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixRHEOTension::init_list(int /*id*/, NeighList *ptr)
|
||||
{
|
||||
list = ptr;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixRHEOTension::setup(int vflag)
|
||||
{
|
||||
// Grow and populate arrays
|
||||
post_force(vflag);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Calculate and apply tension forces
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixRHEOTension::post_force(int vflag)
|
||||
{
|
||||
int i, j, a, ii, jj, inum, jnum, itype, jtype;
|
||||
int fluidi, fluidj;
|
||||
double xtmp, ytmp, ztmp, w, wp, c;
|
||||
double rhoi, rhoj, Voli, Volj;
|
||||
double *dWij, *dWji;
|
||||
double dx[3], ft[3];
|
||||
|
||||
int *ilist, *jlist, *numneigh, **firstneigh;
|
||||
double imass, jmass, rsq, r, rinv;
|
||||
|
||||
int nlocal = atom->nlocal;
|
||||
int newton = force->newton;
|
||||
int dim = domain->dimension;
|
||||
|
||||
v_init(vflag);
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double *rho = atom->rho;
|
||||
double *mass = atom->mass;
|
||||
imageint *image = atom->image;
|
||||
int *type = atom->type;
|
||||
int *status = atom->status;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
if (nmax_store <= atom->nmax)
|
||||
grow_arrays(atom->nmax);
|
||||
|
||||
for (i = 0; i < nlocal+atom->nghost; i++) {
|
||||
cgradt[i][0] = 0.0;
|
||||
cgradt[i][1] = 0.0;
|
||||
cgradt[i][2] = 0.0;
|
||||
norm[i] = 0.0;
|
||||
divnt[i] = 0.0;
|
||||
}
|
||||
|
||||
// Calculate color gradient
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
itype = type[i];
|
||||
fluidi = !(status[i] & PHASECHECK);
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
imass = mass[itype];
|
||||
|
||||
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 = lensq3(dx);
|
||||
|
||||
if (rsq > hsq) continue;
|
||||
|
||||
fluidj = !(status[j] & PHASECHECK);
|
||||
jtype = type[j];
|
||||
r = sqrt(rsq);
|
||||
rinv = 1 / r;
|
||||
|
||||
rhoi = rho[i];
|
||||
rhoj = rho[j];
|
||||
|
||||
// Add corrections for walls
|
||||
if (interface_flag) {
|
||||
if (fluidi && (!fluidj)) {
|
||||
rhoj = compute_interface->correct_rho(j, i);
|
||||
} else if ((!fluidi) && fluidj) {
|
||||
rhoi = compute_interface->correct_rho(i, j);
|
||||
} else if ((!fluidi) && (!fluidj)) {
|
||||
rhoi = rho0;
|
||||
rhoj = rho0;
|
||||
}
|
||||
}
|
||||
|
||||
Voli = mass[itype] / rhoi;
|
||||
Volj = mass[jtype] / rhoj;
|
||||
|
||||
wp = compute_kernel->calc_dw(i, j, dx[0], dx[1], dx[2],r);
|
||||
dWij = compute_kernel->dWij;
|
||||
dWji = compute_kernel->dWji;
|
||||
|
||||
c = 0;
|
||||
if (itype == jtype) c += rhoi;
|
||||
c /= (rhoi + rhoj);
|
||||
|
||||
for (a = 0; a < 3; a++) {
|
||||
cgradt[i][a] -= c * Volj * dWij[a];
|
||||
if (newton || j < nlocal)
|
||||
cgradt[j][a] -= c * Voli * dWji[a];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
comm_stage = 0;
|
||||
comm_reverse = 3;
|
||||
if (newton) comm->reverse_comm(this);
|
||||
|
||||
// Calculate normal direction
|
||||
double minv;
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
minv = 1.0 / sqrt(cgradt[i][0] * cgradt[i][0] + cgradt[i][1] * cgradt[i][1] + cgradt[i][2] * cgradt[i][2]);
|
||||
|
||||
for (a = 0; a < 3; a++)
|
||||
nt[i][a] = cgradt[i][a] * minv;
|
||||
}
|
||||
|
||||
comm->forward_comm(this);
|
||||
|
||||
// Calculate divergence
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
itype = type[i];
|
||||
fluidi = !(status[i] & PHASECHECK);
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
imass = mass[itype];
|
||||
|
||||
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 = lensq3(dx);
|
||||
|
||||
if (rsq > hsq) continue;
|
||||
|
||||
fluidj = !(status[j] & PHASECHECK);
|
||||
jtype = type[j];
|
||||
r = sqrt(rsq);
|
||||
rinv = 1 / r;
|
||||
|
||||
rhoi = rho[i];
|
||||
rhoj = rho[j];
|
||||
|
||||
// Add corrections for walls
|
||||
if (interface_flag) {
|
||||
if (fluidi && (!fluidj)) {
|
||||
rhoj = compute_interface->correct_rho(j, i);
|
||||
} else if ((!fluidi) && fluidj) {
|
||||
rhoi = compute_interface->correct_rho(i, j);
|
||||
} else if ((!fluidi) && (!fluidj)) {
|
||||
rhoi = rho0;
|
||||
rhoj = rho0;
|
||||
}
|
||||
}
|
||||
|
||||
Voli = mass[itype] / rhoi;
|
||||
Volj = mass[jtype] / rhoj;
|
||||
|
||||
wp = compute_kernel->calc_dw(i, j, dx[0], dx[1], dx[2],r);
|
||||
dWij = compute_kernel->dWij;
|
||||
dWji = compute_kernel->dWji;
|
||||
|
||||
for (a = 0; a < 3; a++) {
|
||||
divnt[i] -= nt[i][a] * Volj * dWij[a];
|
||||
norm[i] -= dx[a] * Volj * dWij[a];
|
||||
if (newton || j < nlocal) {
|
||||
divnt[j] -= nt[j][a] * Voli * dWji[a];
|
||||
norm[j] += dx[a] * Voli * dWji[a];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
comm_stage = 1;
|
||||
comm_reverse = 2;
|
||||
if (newton) comm->reverse_comm(this);
|
||||
|
||||
// Skip forces if it's setup
|
||||
if (update->setupflag) return;
|
||||
|
||||
// apply force
|
||||
int prefactor;
|
||||
double unwrap[3];
|
||||
double v[6];
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
itype = type[i];
|
||||
divnt[i] /= norm[i];
|
||||
|
||||
prefactor *= -alpha * divnt[i] / mass[itype];
|
||||
|
||||
for (a = 0; a < 3; a++)
|
||||
f[i][a] += prefactor * cgradt[i][a];
|
||||
|
||||
if (evflag) {
|
||||
domain->unmap(x[i], image[i], unwrap);
|
||||
v[0] = prefactor * cgradt[i][0] * unwrap[0];
|
||||
v[1] = prefactor * cgradt[i][1] * unwrap[1];
|
||||
v[2] = prefactor * cgradt[i][2] * unwrap[2];
|
||||
v[3] = prefactor * cgradt[i][0] * unwrap[1];
|
||||
v[4] = prefactor * cgradt[i][0] * unwrap[2];
|
||||
v[5] = prefactor * cgradt[i][1] * unwrap[2];
|
||||
v_tally(i, v);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
if (evflag) {
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixRHEOTension::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/)
|
||||
{
|
||||
int i, j, a, m;
|
||||
m = 0;
|
||||
|
||||
for (i = 0; i < n; i++) {
|
||||
j = list[i];
|
||||
for (a = 0; a < 3; a++)
|
||||
buf[m++] = nt[j][a];
|
||||
}
|
||||
return m;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixRHEOTension::unpack_forward_comm(int n, int first, double *buf)
|
||||
{
|
||||
int i, a, m, last;
|
||||
|
||||
m = 0;
|
||||
last = first + n;
|
||||
for (i = first; i < last; i++)
|
||||
for (a = 0; a < 3; a++)
|
||||
nt[i][a] = buf[m++];
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixRHEOTension::pack_reverse_comm(int n, int first, double *buf)
|
||||
{
|
||||
int i, a, m, last;
|
||||
|
||||
m = 0;
|
||||
last = first + n;
|
||||
if (comm_stage == 0)
|
||||
for (i = first; i < last; i++)
|
||||
for (a = 0; a < 3; a++)
|
||||
buf[m++] = cgradt[i][a];
|
||||
else
|
||||
for (i = first; i < last; i++) {
|
||||
buf[m++] = norm[i];
|
||||
buf[m++] = divnt[i];
|
||||
}
|
||||
return m;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixRHEOTension::unpack_reverse_comm(int n, int *list, double *buf)
|
||||
{
|
||||
int i, j, a, m;
|
||||
|
||||
m = 0;
|
||||
if (comm_stage == 0)
|
||||
for (i = 0; i < n; i++) {
|
||||
j = list[i];
|
||||
for (a = 0; a < 3; a++)
|
||||
cgradt[j][a] += buf[m++];
|
||||
}
|
||||
else
|
||||
for (i = 0; i < n; i++) {
|
||||
j = list[i];
|
||||
norm[j] += buf[m++];
|
||||
divnt[j] += buf[m++];
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixRHEOTension::grow_arrays(int nmax)
|
||||
{
|
||||
// Grow atom variables and reassign pointers
|
||||
memory->grow(atom->darray[index_cgradt], nmax, 3, "atom:rheo_cgradt");
|
||||
memory->grow(atom->darray[index_nt], nmax, 3, "atom:rheo_nt");
|
||||
memory->grow(atom->dvector[index_divnt], nmax, "atom:rheo_divnt");
|
||||
|
||||
cgradt = atom->darray[index_cgradt];
|
||||
nt = atom->darray[index_nt];
|
||||
divnt = atom->dvector[index_divnt];
|
||||
|
||||
// Grow local variables
|
||||
memory->grow(norm, nmax, "rheo/tension:norm");
|
||||
|
||||
nmax_store = atom->nmax;
|
||||
}
|
||||
@ -11,44 +11,45 @@
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef PAIR_CLASS
|
||||
#ifdef FIX_CLASS
|
||||
// clang-format off
|
||||
PairStyle(rheo/tension,PairRHEOTension)
|
||||
FixStyle(rheo/tension,FixRHEOTension)
|
||||
// clang-format on
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_RHEO_TENSION_H
|
||||
#define LMP_PAIR_RHEO_TENSION_H
|
||||
#ifndef LMP_FIX_RHEO_TENSION_H
|
||||
#define LMP_FIX_RHEO_TENSION_H
|
||||
|
||||
#include "pair.h"
|
||||
#include "fix.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairRHEOTension : public Pair {
|
||||
class FixRHEOTension : public Fix {
|
||||
public:
|
||||
PairRHEOTension(class LAMMPS *);
|
||||
~PairRHEOTension() override;
|
||||
void compute(int, int) override;
|
||||
void settings(int, char **) override;
|
||||
void coeff(int, char **) override;
|
||||
void setup() override;
|
||||
void init_style() override;
|
||||
double init_one(int, int) override;
|
||||
FixRHEOTension(class LAMMPS *, int, char **);
|
||||
~FixRHEOTension() override;
|
||||
int setmask() override;
|
||||
void init() override;
|
||||
void init_list(int, class NeighList *) override;
|
||||
void setup(int) override;
|
||||
void post_force(int) 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;
|
||||
void grow_arrays(int) override;
|
||||
|
||||
protected:
|
||||
int nmax_store;
|
||||
double **nt, *ct;
|
||||
double **alpha;
|
||||
double h, hsq, hinv, hinv3;
|
||||
private:
|
||||
int nmax_store, comm_stage, interface_flag;
|
||||
int index_nt, index_cgradt, index_divnt;
|
||||
|
||||
void allocate();
|
||||
double **nt, **cgradt, *divnt, *norm;
|
||||
double alpha, h, hsq, hinv, hinv3, rho0;
|
||||
|
||||
class ComputeRHEOKernel *compute_kernel;
|
||||
class ComputeRHEOInterface *compute_interface;
|
||||
class FixRHEO *fix_rheo;
|
||||
class NeighList *list;
|
||||
};
|
||||
|
||||
} // namespace LAMMPS_NS
|
||||
@ -1,383 +0,0 @@
|
||||
/* ----------------------------------------------------------------------
|
||||
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)
|
||||
----------------------------------------------------------------------- */
|
||||
|
||||
#include "pair_rheo_tension.h"
|
||||
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "compute_rheo_kernel.h"
|
||||
#include "domain.h"
|
||||
#include "error.h"
|
||||
#include "fix_rheo.h"
|
||||
#include "force.h"
|
||||
#include "math_extra.h"
|
||||
#include "memory.h"
|
||||
#include "modify.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "update.h"
|
||||
#include "utils.h"
|
||||
|
||||
#include <cmath>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace RHEO_NS;
|
||||
using namespace MathExtra;
|
||||
|
||||
static constexpr double EPSILON = 1e-2;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairRHEOTension::PairRHEOTension(LAMMPS *lmp) :
|
||||
Pair(lmp), compute_kernel(nullptr), fix_rheo(nullptr)
|
||||
{
|
||||
restartinfo = 0;
|
||||
single_enable = 0;
|
||||
|
||||
comm_forward = 3;
|
||||
comm_reverse = 3;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairRHEOTension::~PairRHEOTension()
|
||||
{
|
||||
// Remove custom property if it exists
|
||||
int tmp1, tmp2, index;
|
||||
|
||||
index = atom->find_custom("rheo_c_tension", tmp1, tmp2);
|
||||
if (index != -1) atom->remove_custom(index, 1, 0);
|
||||
|
||||
index = atom->find_custom("rheo_n_tension", tmp1, tmp2);
|
||||
if (index != -1) atom->remove_custom(index, 1, 3);
|
||||
|
||||
if (allocated) {
|
||||
memory->destroy(alpha);
|
||||
memory->destroy(setflag);
|
||||
memory->destroy(cutsq);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairRHEOTension::compute(int eflag, int vflag)
|
||||
{
|
||||
int i, j, a, b, ii, jj, inum, jnum, itype, jtype;
|
||||
int fluidi, fluidj;
|
||||
double xtmp, ytmp, ztmp, w, wp;
|
||||
double rhoi, rhoj, voli, volj;
|
||||
double *dWij, *dWji;
|
||||
double dx[3], ft[3];
|
||||
|
||||
int *ilist, *jlist, *numneigh, **firstneigh;
|
||||
double imass, jmass, rsq, r, rinv;
|
||||
|
||||
int nlocal = atom->nlocal;
|
||||
int newton_pair = force->newton_pair;
|
||||
int dim = domain->dimension;
|
||||
|
||||
ev_init(eflag, vflag);
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double *rho = atom->rho;
|
||||
double *mass = atom->mass;
|
||||
double *special_lj = force->special_lj;
|
||||
int *type = atom->type;
|
||||
int *status = atom->status;
|
||||
tagint *tag = atom->tag;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
/*
|
||||
int nmax = atom->nmax;
|
||||
if (nmax_store <= nmax) {
|
||||
memory->grow(ct, nmax, "atom:rheo_c_tension");
|
||||
memory->grow(nnt_tension, nmax, 3, "atom:rheo_n_tension");
|
||||
nmax_store = atom->nmax;
|
||||
}
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
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];
|
||||
imass = mass[itype];
|
||||
rhoi = rho[i];
|
||||
voli = imass / rhoi;
|
||||
fluidi = !(status[i] & PHASECHECK);
|
||||
|
||||
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 = lensq3(dx);
|
||||
jtype = type[j];
|
||||
|
||||
if (rsq > hsq) continue;
|
||||
|
||||
r = sqrt(rsq);
|
||||
rinv = 1 / r;
|
||||
|
||||
jmass = mass[jtype];
|
||||
rhoj = rho[j];
|
||||
volj = jmass / rhoj;
|
||||
fluidj = !(status[j] & PHASECHECK);
|
||||
|
||||
wp = compute_kernel->calc_dw(i, j, dx[0], dx[1], dx[2],r);
|
||||
dWij = compute_kernel->dWij;
|
||||
dWji = compute_kernel->dWji;
|
||||
|
||||
f[i][0] += ft[0];
|
||||
f[i][1] += ft[1];
|
||||
f[i][2] += ft[2];
|
||||
|
||||
if (evflag) // Does not account for unbalanced forces
|
||||
ev_tally_xyz(i, j, nlocal, newton_pair, 0.0, 0.0, ft[0], ft[1], ft[2], dx[0], dx[1], dx[2]);
|
||||
|
||||
if (newton_pair || j < nlocal) {
|
||||
|
||||
f[j][0] -= ft[0];
|
||||
f[j][1] -= ft[1];
|
||||
f[j][2] -= ft[2];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (vflag_fdotr) virial_fdotr_compute();
|
||||
|
||||
comm->reverse_comm(this);
|
||||
comm->forward_comm(this);
|
||||
*/
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
allocate all arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairRHEOTension::allocate()
|
||||
{
|
||||
allocated = 1;
|
||||
int n = atom->ntypes;
|
||||
|
||||
memory->create(setflag, n + 1, n + 1, "pair:setflag");
|
||||
for (int i = 1; i <= n; i++)
|
||||
for (int j = i; j <= n; j++)
|
||||
setflag[i][j] = 0;
|
||||
|
||||
memory->create(alpha, n + 1, n + 1, "pair:alpha");
|
||||
memory->create(cutsq, n + 1, n + 1, "pair:cutsq");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairRHEOTension::settings(int narg, char **arg)
|
||||
{
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more type pairs
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairRHEOTension::coeff(int narg, char **arg)
|
||||
{
|
||||
if (narg != 3)
|
||||
error->all(FLERR,"Incorrect number of args for pair_style rheo coefficients");
|
||||
if (!allocated)
|
||||
allocate();
|
||||
|
||||
int ilo, ihi, jlo, jhi;
|
||||
utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi,error);
|
||||
utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi,error);
|
||||
|
||||
double alpha_one = utils::numeric(FLERR, arg[2], false, lmp);
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = 0; j <= atom->ntypes; j++) {
|
||||
alpha[i][j] = alpha_one;
|
||||
setflag[i][j] = 1;
|
||||
count++;
|
||||
}
|
||||
}
|
||||
|
||||
if (count == 0)
|
||||
error->all(FLERR,"Incorrect args for pair rheo/tension coefficients");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
setup specific to this pair style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairRHEOTension::setup()
|
||||
{
|
||||
auto fixes = modify->get_fix_by_style("rheo");
|
||||
if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use pair rheo");
|
||||
fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
|
||||
/*
|
||||
compute_kernel = fix_rheo->compute_kernel;
|
||||
compute_grad = fix_rheo->compute_grad;
|
||||
compute_interface = fix_rheo->compute_interface;
|
||||
h = fix_rheo->h;
|
||||
csq = fix_rheo->csq;
|
||||
rho0 = fix_rheo->rho0;
|
||||
|
||||
hsq = h * h;
|
||||
hinv = 1.0 / h;
|
||||
hinv3 = hinv * 3.0;
|
||||
*/
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init specific to this pair style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairRHEOTension::init_style()
|
||||
{
|
||||
neighbor->add_request(this);
|
||||
|
||||
|
||||
// Create c_tension arrays n_tension arrays if they don'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
|
||||
// For B and gradC, create a local array since they are unlikely to be printed
|
||||
|
||||
int tmp1, tmp2;
|
||||
int index = atom->find_custom("rheo_c_tension", tmp1, tmp2);
|
||||
if (index == -1) index = atom->add_custom("rheo_c_tension", 1, 0);
|
||||
ct = atom->dvector[index];
|
||||
|
||||
index = atom->find_custom("rheo_n_tension", tmp1, tmp2);
|
||||
if (index == -1) index = atom->add_custom("rheo_n_tension", 1, 3);
|
||||
nt = atom->darray[index];
|
||||
|
||||
nmax_store = atom->nmax;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init for one type pair i,j and corresponding j,i
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairRHEOTension::init_one(int i, int j)
|
||||
{
|
||||
if (setflag[i][j] == 0)
|
||||
error->all(FLERR,"All pair rheo/tension coeffs are not set");
|
||||
|
||||
alpha[j][i] = alpha[i][j];
|
||||
|
||||
return h;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int PairRHEOTension::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/)
|
||||
{
|
||||
/*
|
||||
int i,j,k,m;
|
||||
m = 0;
|
||||
double *rho = atom->rho;
|
||||
|
||||
for (i = 0; i < n; i++) {
|
||||
j = list[i];
|
||||
if (comm_stage == 0) {
|
||||
buf[m++] = fp_store[j][0];
|
||||
buf[m++] = fp_store[j][1];
|
||||
buf[m++] = fp_store[j][2];
|
||||
} else {
|
||||
buf[m++] = chi[j];
|
||||
buf[m++] = rho[j];
|
||||
}
|
||||
}
|
||||
return m;
|
||||
*/
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairRHEOTension::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++) {
|
||||
if (comm_stage == 0) {
|
||||
fp_store[i][0] = buf[m++];
|
||||
fp_store[i][1] = buf[m++];
|
||||
fp_store[i][2] = buf[m++];
|
||||
} else {
|
||||
chi[i] = buf[m++];
|
||||
rho[i] = buf[m++];
|
||||
}
|
||||
}
|
||||
*/
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int PairRHEOTension::pack_reverse_comm(int n, int first, double *buf)
|
||||
{
|
||||
/*
|
||||
int i, k, m, last;
|
||||
double **fp_store = compute_interface->fp_store;
|
||||
|
||||
m = 0;
|
||||
last = first + n;
|
||||
for (i = first; i < last; i++) {
|
||||
buf[m++] = fp_store[i][0];
|
||||
buf[m++] = fp_store[i][1];
|
||||
buf[m++] = fp_store[i][2];
|
||||
}
|
||||
|
||||
return m;
|
||||
*/
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairRHEOTension::unpack_reverse_comm(int n, int *list, double *buf)
|
||||
{
|
||||
/*
|
||||
int i, j, k, m;
|
||||
double **fp_store = compute_interface->fp_store;
|
||||
|
||||
m = 0;
|
||||
for (i = 0; i < n; i++) {
|
||||
j = list[i];
|
||||
fp_store[j][0] += buf[m++];
|
||||
fp_store[j][1] += buf[m++];
|
||||
fp_store[j][2] += buf[m++];
|
||||
}
|
||||
*/
|
||||
}
|
||||
Reference in New Issue
Block a user