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

Initial surface compute

Browse Source
This commit is contained in:
jtclemm
2023-04-19 17:15:00 -06:00
parent 5980fdf9fd
commit d85ce6a392
9 changed files with 714 additions and 25 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/RHEO/compute_rheo_grad.h
View File

@ -13,7 +13,7 @@
#ifdef COMPUTE_CLASS
// clang-format off
ComputeStyle(rheo/grad,ComputeRHEOGrad)
ComputeStyle(RHEO/GRAD,ComputeRHEOGrad)
// clang-format on
#else

2
src/RHEO/compute_rheo_interface.cpp
View File

@ -398,7 +398,7 @@ void ComputeRHEOInterface::store_forces()
double ComputeRHEOInterface::memory_usage()
{
double bytes = nmax * sizeof(double);
double bytes = 3 * nmax_old * sizeof(double);
return bytes;
}

2
src/RHEO/compute_rheo_interface.h
View File

@ -13,7 +13,7 @@
#ifdef COMPUTE_CLASS
// clang-format off
ComputeStyle(rheo/interface,ComputeRHEOInterface)
ComputeStyle(RHEO/INTERFACE,ComputeRHEOInterface)
// clang-format on
#else

2
src/RHEO/compute_rheo_kernel.h
View File

@ -13,7 +13,7 @@
#ifdef COMPUTE_CLASS
// clang-format off
ComputeStyle(rheo/kernel,ComputeRHEOKernel)
ComputeStyle(RHEO/KERNEL,ComputeRHEOKernel)
// clang-format on
#else

595
src/RHEO/compute_rheo_surface.cpp Normal file
View File

@ -0,0 +1,595 @@
/* ----------------------------------------------------------------------
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_surface.h"
#include "fix_rheo.h"
#include "compute_rheo_kernel.h"
#include "compute_rheo_solids.h"
#include "atom.h"
#include "memory.h"
#include "atom.h"
#include "comm.h"
#include "modify.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "error.h"
#include "force.h"
#include "domain.h"
using namespace LAMMPS_NS;
using namespace FixConst;
/* ---------------------------------------------------------------------- */
ComputeRHEOSurface::ComputeRHEOSurface(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 6) error->all(FLERR,"Illegal fix rheo/surface command");
cut = utils::numeric(FLERR,arg[3],false,lmp);
divR_limit = utils::numeric(FLERR,arg[4],false,lmp);
coord_limit = utils::inumeric(FLERR,arg[5],false,lmp);
divr_flag = 1;
if (narg == 7) {
divr_flag = 0;
}
int dim = domain->dimension;
peratom_flag = 1;
size_peratom_cols = dim;
peratom_freq = 1;
comm_forward = 2;
comm_reverse = dim*dim + 1;
cutsq = cut*cut;
B = nullptr;
gradC = nullptr;
n_surface = nullptr;
int nall = atom->nlocal + atom->nghost;
nmax = nall;
memory->create(B,nmax,dim*dim,"fix/rheo/surface:B");
memory->create(gradC,nmax,dim*dim,"fix/rheo/surface:gradC");
memory->create(n_surface,nmax,dim,"fix/rheo/surface:B");
array_atom = n_surface;
compute_kernel = nullptr;
compute_solids = NULL;
fix_rheo = nullptr;
}
/* ---------------------------------------------------------------------- */
ComputeRHEOSurface::~ComputeRHEOSurface()
{
if (modify->nfix) modify->delete_fix("PROPERTY_ATOM_RHEO_SURFACE");
memory->destroy(B);
memory->destroy(gradC);
memory->destroy(n_surface);
}
void ComputeRHEOSurface::post_constructor()
{
//Store persistent per atom quantities
char **fixarg = new char*[5];
fixarg[0] = (char *) "PROPERTY_ATOM_RHEO_SURFACE";
fixarg[1] = (char *) "all";
fixarg[2] = (char *) "property/atom";
fixarg[3] = (char *) "d_divr";
fixarg[4] = (char *) "d_rsurf";
modify->add_fix(5,fixarg,1);
int temp_flag;
index_divr = atom->find_custom("divr", temp_flag);
if ((index_divr < 0) || (temp_flag != 1))
error->all(FLERR, "Pair rheo/surface can't find fix property/atom divr");
index_rsurf = atom->find_custom("rsurf", temp_flag);
if ((index_rsurf < 0) || (temp_flag != 1))
error->all(FLERR, "Pair rheo/surface can't find fix property/atom rsurf");
delete [] fixarg;
divr = atom->dvector[index_divr];
}
/* ---------------------------------------------------------------------- */
int ComputeRHEOSurface::setmask()
{
int mask = 0;
mask |= PRE_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOSurface::init()
{
// need an occasional full neighbor list
int irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->fix = 1;
neighbor->requests[irequest]->half = 1;
neighbor->requests[irequest]->full = 0;
int flag;
int ifix = modify->find_fix_by_style("rheo");
if (ifix == -1) error->all(FLERR, "Need to define fix rheo to use fix rheo/surface");
fix_rheo = ((FixRHEO *) modify->fix[ifix]);
compute_kernel = fix_rheo->compute_kernel;
compute_solids = fix_rheo->compute_solids;
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOSurface::setup_pre_force(int /*vflag*/)
{
pre_force(0);
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOSurface::init_list(int /*id*/, NeighList *ptr)
{
list = ptr;
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOSurface::pre_force(int /*vflag*/)
{
int i, j, ii, jj, jnum, a, b, itype, jtype;
double xtmp, ytmp, ztmp, delx, dely, delz, rsq, r, wp, Voli, Volj, rhoi, rhoj;
int *jlist;
double *dWij, *dWji;
double dx[3];
divr = atom->dvector[index_divr];
// neighbor list variables
int inum, *ilist, *numneigh, **firstneigh;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
double **x = atom->x;
int *surface = atom->surface;
int *phase = atom->phase;
double *rsurf = atom->dvector[index_rsurf];
int newton = force->newton;
int dim = domain->dimension;
int *mask = atom->mask;
int *type = atom->type;
double *mass = atom->mass;
double *rho = atom->rho;
double *temp = atom->temp;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
if (nmax <= nall) {
nmax = nall;
memory->destroy(B);
memory->destroy(gradC);
memory->destroy(n_surface);
memory->create(B,nmax,dim*dim,"fix/rheo/surface:B");
memory->create(gradC,nmax,dim*dim,"fix/rheo/surface:gradC");
memory->create(n_surface,nmax,dim,"fix/rheo/surface:n_surface");
array_atom = n_surface;
}
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;
}
n_surface[i][a] = 0.0;
}
divr[i] = 0.0;
surface[i] = 0;
}
// loop over neighbors to calculate the average orientation of neighbors
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];
itype = type[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];
dx[0] = delx;
dx[1] = dely;
dx[2] = delz;
rsq = delx * delx + dely * dely + delz * delz;
if (rsq < cutsq) {
jtype = type[j];
rhoi = rho[i];
rhoj = rho[j];
// Add corrections for walls
if (phase[i] <= FixRHEO::FLUID_MAX && phase[j] > FixRHEO::FLUID_MAX) {
rhoj = compute_solids->correct_rho(j,i);
} else if (phase[i] > FixRHEO::FLUID_MAX && phase[j] <= FixRHEO::FLUID_MAX) {
rhoi = compute_solids->correct_rho(i,j);
} else if (phase[i] > FixRHEO::FLUID_MAX && phase[j] > FixRHEO::FLUID_MAX) {
rhoi = 1.0;
rhoj = 1.0;
}
Voli = mass[itype]/rhoi;
Volj = mass[jtype]/rhoj;
//compute kernel gradient
wp = compute_kernel->calc_dw_quintic(i, j, delx, dely, delz, sqrt(rsq),compute_kernel->dWij,compute_kernel->dWji);
//wp = compute_kernel->calc_dw(i, j, delx, dely, delz, sqrt(rsq));//,compute_kernel->dWij,compute_kernel->dWji);
dWij = compute_kernel->dWij;
dWji = compute_kernel->dWji;
for (a=0; a<dim; a++){
divr[i] -= dWij[a]*dx[a]*Volj; // dx = xi-xj = xji = -xij
gradC[i][a] += dWij[a]*Volj;
}
if (j < nlocal || newton) {
for (a=0; a<dim; a++){
divr[j] += dWji[a]*dx[a]*Voli;
gradC[j][a] += dWji[a]*Voli;
}
}
}
}
}
comm_stage = 0;
comm_reverse = dim*dim + 1; // gradC and divr
comm_forward = 1; // divr
if (newton) comm->reverse_comm_fix(this);
comm->forward_comm_fix(this);
int *coordination = compute_kernel->coordination;
// Find the free-surface
//0-bulk 1-surf vicinity 2-surface 3-splash
if (divr_flag) {
for (i = 0; i < nall; i++) {
if (mask[i] & groupbit) {
surface[i] = 0;
rsurf[i] = cut; //Maximum range that can be seen
if (divr[i] < divR_limit) {
surface[i] = 2;
rsurf[i] = 0.0;
if (coordination[i] < coord_limit) surface[i] = 3;
}
}
}
} else {
for (i = 0; i < nall; i++) {
if (mask[i] & groupbit) {
surface[i] = 0;
rsurf[i] = cut; //Maximum range that can be seen
if (coordination[i] < divR_limit) {
surface[i] = 2;
rsurf[i] = 0.0;
if (coordination[i] < coord_limit) surface[i] = 3;
}
}
}
}
//comm_stage = 1;
//comm_forward = 1;
//comm->forward_comm_fix(this); // communicate free surface particles
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];
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) {
r = sqrt(rsq);
if (surface[i] == 0 && surface[j] == 2) surface[i] = 1;
if (surface[j] == 0 && surface[i] == 2) surface[j] = 1;
if (surface[j] == 2) rsurf[i] = MIN(rsurf[i], r);
if (surface[i] == 2) rsurf[j] = MIN(rsurf[j], r);
}
}
}
comm_stage = 1;
comm_reverse = 2;
comm_forward = 2;
if (newton) comm->reverse_comm_fix(this);
comm->forward_comm_fix(this);
//Now loop again and for each surface particle (2)
// find its neighbors that are bulk (0) and convert to surface vicinity (1)
// if the surface particle has no (0) or (1) neighbors then it is a spash (3)
//for (ii = 0; ii < inum; ii++) { // is this the right i and j loop for this?
// i = ilist[ii];
//
// if (surface[i]!=2) continue; //Only consider surface particles
//
// bool nobulkneigh = true; // whether we have no bulk neighbors
// xtmp = x[i][0];
// ytmp = x[i][1];
// ztmp = x[i][2];
// jlist = firstneigh[i];
// jnum = numneigh[i];
//
// for (jj = 0; jj < jnum; jj++) {
// j = jlist[jj];
// j &= NEIGHMASK;
//
// //other surface or splash neighbors do not need labeling
// if (surface[j]>=2){
// continue;
// }
//
// //check distance criterion rij < h = cutsq/9 for quintic kernel
// delx = xtmp - x[j][0];
// dely = ytmp - x[j][1];
// delz = ztmp - x[j][2];
// dx[0] = 3.0*delx; // multiplied by three here to make criterion r<h instead of r<3*h
// dx[1] = 3.0*dely;
// dx[2] = 3.0*delz;
// rsq = delx * delx + dely * dely + delz * delz;
// if (rsq < cutsq) {
// //We have identified 1 bulk fluid neighbor
// nobulkneigh = false;
// //that bulk fluid neighbor is in the vicinity of hte surface
// surface[j] = 1;
// }
// }
// if (nobulkneigh){
// surface[i] = 3;
// }
//}
//
// //Reverse comm surface?
//
// // loop over neighbors to calculate the average orientation
// // skip for bulk or splash
// for (ii = 0; ii < inum; ii++) {
// i = ilist[ii];
// if ((surface[i]==0)||(surface[i]==3)){
// continue;
// }
//
// itype = type[i];
// rhoi = rho[i];
// Voli = mass[itype]/rhoi;
//
// xtmp = x[i][0];
// ytmp = x[i][1];
// ztmp = x[i][2];
//
// 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];
// dx[0] = delx; // multiplied by three here to make criterion r<h instead of r<3*h
// dx[1] = dely;
// dx[2] = delz;
// rsq = delx * delx + dely * dely + delz * delz;
// if (rsq < cutsq) {
//
// jtype = type[j];
// rhoj = rho[j];
// Volj = mass[jtype]/rhoj;
//
// for (a=0; a<dim; a++){
// for (b=0; b<dim; b++){
// B[i][a*dim+b] -= dx[a]*dWij[b]*Volj;
// }
// }
//
// if (j < nlocal || newton) {
// for (a=0; a<dim; a++){
// for (b=0; b<dim; b++){
// B[j][a*dim+b] += dx[a]*dWji[b]*Voli;
// }
// }
// }
// }
// }
// }
// //reverse comm to populate B[j] if Newton is on
// comm_stage = 2;
// comm_reverse = dim*dim; // B
// if (newton) comm->reverse_comm_fix(this);
//
//
// // Now need to invert each B
// int status, s;
// //LU requires a permuation matrix
// gsl_permutation * p = gsl_permutation_alloc(dim);
// for (ii = 0; ii < inum; ii++) {
// i = ilist[ii];
// if ((surface[i]==0)||(surface[i]==3)){
// continue;
// }
//
// //Use gsl to get Binv
// //B is not symmteric so we will use a LU decomp
// gsl_matrix_view gB = gsl_matrix_view_array(B[i],dim,dim);
// status = 0;
// status = gsl_linalg_LU_decomp(&gB.matrix,p,&s); //B[i] is now the LU decomp
// // check if decomposition failure
// if (status) {
// fprintf(stderr, "failed, gsl_errno=%d.n", status);
// continue;
// } else {
// gsl_linalg_LU_invx(&gB.matrix,p); //B[i] is now inv(B[i])
// }
// }
// gsl_permutation_free(p);
double maggC = 0.0;
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
maggC=0;
for (a=0;a<dim;a++){
maggC += gradC[i][a]*gradC[i][a];
}
maggC = sqrt(maggC) + 1e-10;
for (a=0;a<dim;a++){
n_surface[i][a] = -gradC[i][a]/maggC;
}//dr can then be calculated by fix vshift
}
}
}
/* ---------------------------------------------------------------------- */
int ComputeRHEOSurface::pack_reverse_comm(int n, int first, double *buf)
{
int i,a,b,k,m,last;
int dim = domain->dimension;
int *surface = atom->surface;
double *rsurf = atom->dvector[index_rsurf];
m = 0;
last = first + n;
for (i = first; i < last; i++) {
if (comm_stage == 0) {
buf[m++] = divr[i];
for (a = 0; a < dim; a ++ )
for (b = 0; b < dim; b ++)
buf[m++] = gradC[i][a*dim + b];
} else if (comm_stage == 1) {
buf[m++] = (double) surface[i];
buf[m++] = rsurf[i];
} else if (comm_stage == 2) {
for (a = 0; a < dim; a ++ )
for (b = 0; b < dim; b ++)
buf[m++] = B[i][a*dim + b];
}
}
return m;
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOSurface::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,a,b,k,j,m;
int dim = domain->dimension;
int *surface = atom->surface;
double *rsurf = atom->dvector[index_rsurf];
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
if (comm_stage == 0) {
divr[j] += buf[m++];
for (a = 0; a < dim; a ++ )
for (b = 0; b < dim; b ++)
gradC[j][a*dim + b] += buf[m++];
} else if (comm_stage == 1) {
int temp = (int) buf[m++];
surface[j] = MAX(surface[j], temp);
double temp2 = buf[m++];
rsurf[j] = MIN(rsurf[j], temp2);
} else if (comm_stage == 2) {
for (a = 0; a < dim; a ++ )
for (b = 0; b < dim; b ++)
B[j][a*dim + b] += buf[m++];
}
}
}
/* ---------------------------------------------------------------------- */
int ComputeRHEOSurface::pack_forward_comm(int n, int *list, double *buf,
int /*pbc_flag*/, int * /*pbc*/)
{
int i,j,a,b,k,m;
int *surface = atom->surface;
double *rsurf = atom->dvector[index_rsurf];
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
if (comm_stage == 0) {
buf[m++] = divr[j];
} else if (comm_stage == 1) {
buf[m++] = (double) surface[j];
buf[m++] = rsurf[j];
}
}
return m;
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOSurface::unpack_forward_comm(int n, int first, double *buf)
{
int i, k, a, b, m, last;
int *surface = atom->surface;
double *rsurf = atom->dvector[index_rsurf];
m = 0;
last = first + n;
for (i = first; i < last; i++) {
if (comm_stage == 0) {
divr[i] = buf[m++];
} else if (comm_stage == 1) {
surface[i] = (int) buf[m++];
rsurf[i] = buf[m++];
}
}
}

72
src/RHEO/compute_rheo_surface.h Normal file
View File

@ -0,0 +1,72 @@
/* -*- 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/SURFACE,ComputeRHEOSurface)
// clang-format on
#else
#ifndef LMP_COMPUTE_RHEO_INTERFACE_H
#define LMP_COMPUTE_RHEO_INTERFACE_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeRHEOSurface : public Compute {
public:
ComputeRHEOSurface(class LAMMPS *, int, char **);
~ComputeRHEOSurface();
void init() override;
void init_list(int, class NeighList *) override;
void compute_peratom() override;
int pack_reverse_comm(int, int, double *) override;
void unpack_reverse_comm(int, int *, double *) override;
int pack_forward_comm(int, int *, double *, int, int *) override;
void unpack_forward_comm(int, int, double *) override;
double **gradC, **n_surface;
private:
double cut, cutsq, threshold;
int surface_style, nmax_old;
double **B, *divr;
int comm_stage;
int index_divr;
int index_rsurf;
double divR_limit;
int coord_limit;
class NeighList *list;
class FixRHEO *fix_rheo;
class ComputeRHEOKernel *compute_kernel;
class ComputeRHEOSolids *compute_solids;
};
}
#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.
*/

2
src/RHEO/compute_rheo_vshift.h
View File

@ -13,7 +13,7 @@
#ifdef COMPUTE_CLASS
// clang-format off
ComputeStyle(rheo/vshift,ComputeRHEOVShift)
ComputeStyle(RHEO/VSHIFT,ComputeRHEOVShift)
// clang-format on
#else

52
src/RHEO/fix_rheo.cpp
View File

@ -88,12 +88,24 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
thermal_flag = 1;
} else if (strcmp(arg[iarg],"surface/detection") == 0) {
surface_flag = 1;
if(iarg + 2 >= narg) error->all(FLERR,"Illegal surface/detection option in fix rheo");
if (strcmp(arg[iarg + 1], "coordination")) {
surface_style = COORDINATION;
zmin_surface = utils::inumeric(FLERR,arg[iarg + 2],false,lmp);
} else if (strcmp(arg[iarg + 1], "divergence")) {
surface_style = DIVR;
divr_surface = utils::numeric(FLERR,arg[iarg + 2],false,lmp);
} else {
error->all(FLERR,"Illegal surface/detection option in fix rheo, {}", arg[iarg + 1]);
}
iarg += 2;
} else if (strcmp(arg[iarg],"interface/reconstruction") == 0) {
interface_flag = 1;
} else if (strcmp(arg[iarg],"rhosum") == 0) {
rhosum_flag = 1;
if(iarg + 1 >= narg) error->all(FLERR,"Illegal rhosum option in fix rheo");
zmin_rhosum = utils::inumeric(FLERR,arg[iarg + 1],false,lmp);
rhosum_zmin = utils::inumeric(FLERR,arg[iarg + 1],false,lmp);
iarg += 1;
} else if (strcmp(arg[iarg],"rho0") == 0) {
if(iarg + 1 >= narg) error->all(FLERR,"Illegal rho0 option in fix rheo");
@ -117,6 +129,7 @@ FixRHEO::~FixRHEO()
if (compute_kernel) modify->delete_compute("rheo_kernel");
if (compute_grad) modify->delete_compute("rheo_grad");
if (compute_interface) modify->delete_compute("rheo_interface");
if (compute_surface) modify->delete_compute("compute_surface");
if (compute_rhosum) modify->delete_compute("rheo_rhosum");
if (compute_vshift) modify->delete_compute("rheo_vshift");
}
@ -128,28 +141,33 @@ 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"));
compute_kernel->fix_rheo = this;
std::string cmd = "rheo_grad all rheo/grad velocity rho viscosity";
std::string cmd = "rheo_grad all RHEO/GRAD velocity rho viscosity";
if (thermal_flag) cmd += "temperature";
compute_grad = dynamic_cast<ComputeRHEOGrad *>(modify->add_compute(cmd));
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_rhosum all RHEO/RHO/SUM"));
compute_rhosum->fix_rheo = this;
}
if (shift_flag) {
compute_vshift = dynamic_cast<ComputeRHEOVShift *>(modify->add_compute("rheo_vshift all rheo/vshift"));
compute_vshift = dynamic_cast<ComputeRHEOVShift *>(modify->add_compute("rheo_vshift all RHEO/VSHIFT"));
compute_vshift->fix_rheo = this;
}
if (interface_flag) {
compute_interface = dynamic_cast<ComputeRHEOInterface *>(modify->add_compute(fmt::format("rheo_interface all rheo/interface")));
compute_interface = dynamic_cast<ComputeRHEOInterface *>(modify->add_compute(fmt::format("rheo_interface all RHEO/INTERFACE")));
compute_interface->fix_rheo = this;
}
if (surface_flag) {
compute_surface = dynamic_cast<ComputeRHEOSurface *>(modify->add_compute(fmt::format("rheo_surface all RHEO/SURFACE")));
compute_surface->fix_rheo = this;
}
}
/* ---------------------------------------------------------------------- */
@ -180,10 +198,11 @@ void FixRHEO::setup_pre_force(int /*vflag*/)
{
// Check to confirm accessory fixes do not preceed FixRHEO
// Note: these fixes set this flag in setup_pre_force()
if (viscosity_fix_defined || pressure_fix_defined || thermal_fix_defined || surface_fix_defined)
if (viscosity_fix_defined || pressure_fix_defined || thermal_fix_defined)
error->all(FLERR, "Fix RHEO must be defined before all other RHEO fixes");
pre_force(0);
// Calculate surfaces
if (surface_flag) compute_surface->compute_peratom();
}
/* ---------------------------------------------------------------------- */
@ -201,14 +220,10 @@ void FixRHEO::setup()
if(!thermal_fix_defined && thermal_flag)
error->all(FLERR, "Missing fix rheo/thermal");
if(!surface_fix_defined && surface_flag)
error->all(FLERR, "Missing fix rheo/surface");
// Reset to zero for next run
thermal_fix_defined = 0;
viscosity_fix_defined = 0;
pressure_fix_defined = 0;
surface_fix_defined = 0;
// Check fixes cover all atoms (doesnt ensure user covers atoms created midrun)
// (pressure is currently required to be group all)
@ -259,7 +274,8 @@ void FixRHEO::initial_integrate(int /*vflag*/)
double **gradr = compute_grad->gradr;
double **gradv = compute_grad->gradv;
double **vshift = compute_vshift->vshift;
double **vshift;
if (shift_flag) compute_vshift->vshift;
int *type = atom->type;
int *mask = atom->mask;
@ -287,8 +303,11 @@ void FixRHEO::initial_integrate(int /*vflag*/)
// Update gradients and interpolate solid properties
compute_grad->forward_fields(); // also forwards v and rho for chi
compute_interface->store_forces(); // Need to save, wiped in exchange
if (interface_flag) {
// Need to save, wiped in exchange
compute_interface->store_forces();
compute_interface->compute_peratom();
}
compute_grad->compute_peratom();
// Position half-step
@ -350,7 +369,7 @@ void FixRHEO::pre_force(int /*vflag*/)
compute_grad->forward_fields(); // also forwards v and rho for chi
compute_kernel->compute_peratom();
compute_interface->compute_peratom();
if (interface_flag) compute_interface->compute_peratom();
compute_grad->compute_peratom();
compute_grad->forward_gradients();
@ -369,6 +388,9 @@ void FixRHEO::pre_force(int /*vflag*/)
status[i] &= ~STATUS_SHIFT;
}
}
// Calculate surfaces, update status
if (surface_flag) compute_surface->compute_peratom();
}
/* ---------------------------------------------------------------------- */

8
src/RHEO/fix_rheo.h
View File

@ -40,9 +40,11 @@ class FixRHEO : public Fix {
// Model parameters
double h, rho0, csq;
int zmin_kernel, rhosum_zmin;
int kernel_style;
int zmin_kernel, zmin_rhosum, zmin_surface;
int kernel_style, surface_style;
double divr_surface;
enum {QUINTIC, CRK0, CRK1, CRK2};
enum {COORDINATION, DIVR};
// Status variables
enum {
@ -75,8 +77,6 @@ class FixRHEO : public Fix {
int viscosity_fix_defined;
int pressure_fix_defined;
int thermal_fix_defined;
int interface_fix_defined;
int surface_fix_defined;
class ComputeRHEOGrad *compute_grad;
class ComputeRHEOKernel *compute_kernel;