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

Updating surface compute

Browse Source
This commit is contained in:
jtclemm
2023-04-20 14:45:35 -06:00
parent d85ce6a392
commit a4d971df52
9 changed files with 195 additions and 373 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/RHEO/compute_rheo_grad.cpp
View File

@ -38,7 +38,7 @@ enum{COMMGRAD, COMMFIELD};
/* ---------------------------------------------------------------------- */
ComputeRHEOGrad::ComputeRHEOGrad(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), fix_rheo(nullptr), compute_interface(nullptr), compute_kernel(nullptr),
Compute(lmp, narg, arg), fix_rheo(nullptr), list(nullptr), compute_interface(nullptr), compute_kernel(nullptr),
gradv(nullptr), gradr(nullptr), gradt(nullptr), gradn(nullptr)
{
if (narg < 4) error->all(FLERR,"Illegal compute rheo/grad command");

1
src/RHEO/compute_rheo_interface.cpp
View File

@ -84,7 +84,6 @@ void ComputeRHEOInterface::init()
// Do not create grow callback as there's no reason to copy/exchange data
// Manually grow if nmax_old exceeded
int create_flag = 0;
int tmp1, tmp2;
int nmax = atom->nmax;
int index = atom->find_custom("rheo_chi", tmp1, tmp2);

12
src/RHEO/compute_rheo_kernel.cpp
View File

@ -55,7 +55,7 @@ Move away from h notation, use cut?
ComputeRHEOKernel::ComputeRHEOKernel(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
C(nullptr), C0(nullptr), coordination(nullptr), compute_interface(nullptr)
list(nullptr), C(nullptr), C0(nullptr), coordination(nullptr), compute_interface(nullptr)
{
if (narg != 3) error->all(FLERR,"Illegal compute rheo/kernel command");
@ -137,18 +137,18 @@ void ComputeRHEOKernel::init()
ncor = 0;
Mdim = 0;
if (kernel_type == CRK0) {
memory->create(C0, nmax, "rheo/kernel:C0");
memory->create(C0, nmax_old, "rheo/kernel:C0");
} else if (kernel_type == CRK1) {
Mdim = 1 + dim;
ncor = 1 + dim;
memory->create(C, nmax, ncor, Mdim, "rheo/kernel:C");
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, ncor, Mdim, "rheo/kernel:C");
memory->create(C, nmax_old, ncor, Mdim, "rheo/kernel:C");
comm_forward = ncor * Mdim;
}
}
@ -491,7 +491,7 @@ void ComputeRHEOKernel::compute_peratom()
gsl_error_flag = 0;
gsl_error_tags.clear();
int i, j, ii, jj, jnum, g, a, b, gsl_error;
int i, j, ii, jj, inum, jnum, g, a, b, gsl_error;
double xtmp, ytmp, ztmp, r, rsq, w, vj;
double dx[3];
gsl_matrix_view gM;
@ -506,7 +506,7 @@ void ComputeRHEOKernel::compute_peratom()
int *status = atom->status;
tagint *tag = atom->tag;
int inum, *ilist, *jlist, *numneigh, **firstneigh;
int *ilist, *jlist, *numneigh, **firstneigh;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;

528
src/RHEO/compute_rheo_surface.cpp
View File

@ -18,135 +18,97 @@
#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 "compute_rheo_kernel.h"
#include "compute_rheo_solids.h"
#include "domain.h"
#include "error.h"
#include "fix_rheo.h"
#include "force.h"
#include "math_extra.h"
#include "memory.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;
using namespace MathExtra;
#define epsilon 1e-10;
/* ---------------------------------------------------------------------- */
ComputeRHEOSurface::ComputeRHEOSurface(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
Fix(lmp, narg, arg), fix_rheo(nullptr), list(nullptr), compute_kernel(nullptr), compute_solids(nullptr),
B(nullptr), gradC(nullptr), nsurface(nullptr), divr(nullptr), rsurface(nullptr)
{
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;
}
if (narg != 3) error->all(FLERR,"Illegal fix RHEO/SURFACE command");
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;
comm_reverse = dim * dim + 1;
}
/* ---------------------------------------------------------------------- */
ComputeRHEOSurface::~ComputeRHEOSurface()
{
if (modify->nfix) modify->delete_fix("PROPERTY_ATOM_RHEO_SURFACE");
// Remove custom property if it exists
int tmp1, tmp2, index;
index = atom->find_custom("rheo_divr", tmp1, tmp2);
if (index != -1) atom->remove_custom(index, 1, 0);
index = atom->find_custom("rheo_rsurface", tmp1, tmp2);
if (index != -1) atom->remove_custom(index, 1, 0);
index = atom->find_custom("rheo_nsurface", tmp1, tmp2);
if (index != -1) atom->remove_custom(index, 1, 3);
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;
}
cut = fix_rheo->cut;
rho0 = fix_rheo->rho0;
threshold_style = fix_rheo->surface_style;
threshold_divr = fix_rheo->divrsurface;
threshold_z = fix_rheo->zminsurface;
/* ---------------------------------------------------------------------- */
cutsq = cut * cut;
void ComputeRHEOSurface::setup_pre_force(int /*vflag*/)
{
pre_force(0);
// 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_old 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_divr", tmp1, tmp2);
if (index == -1) index = atom->add_custom("rheo_divr", 1, 0);
divr = atom->dvector[index];
index = atom->find_custom("rheo_rsurface", tmp1, tmp2);
if (index == -1) index = atom->add_custom("rheo_rsurface", 1, 0);
rsurface = atom->dvector[index];
index = atom->find_custom("rheo_nsurface", tmp1, tmp2);
if (index == -1) index = atom->add_custom("rheo_nsurface", 1, 3);
nsurface = atom->darray[index];
nmax_old = atom->nmax;
memory->create(B, nmax_old, dim * dim, "rheo/surface:B");
memory->create(gradC, nmax_old, dim * dim, "rheo/surface:gradC");
// need an occasional half neighbor list
neighbor->add_request(this, NeighConst::REQ_HALF);
}
/* ---------------------------------------------------------------------- */
@ -158,60 +120,56 @@ void ComputeRHEOSurface::init_list(int /*id*/, NeighList *ptr)
/* ---------------------------------------------------------------------- */
void ComputeRHEOSurface::pre_force(int /*vflag*/)
void ComputeRHEOSurface::compute_peratom()
{
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;
int i, j, ii, jj, inum, jnum, a, b, itype, jtype, fluidi, fluidj;
double xtmp, ytmp, ztmp, rsq, Voli, Volj, rhoi, rhoj;
double *dWij, *dWji;
double dx[3];
int *ilist, *jlist, *numneigh, **firstneigh;
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 *status = atom->status;
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;
int *coordination = compute_kernel->coordination;
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);
int nmax = atom->nmax;
if (nmax_old <= nmax) {
memory->grow(divr, nmax, "atom:rheo_divr");
memory->grow(rsurface, nmax, "atom:rheo_rsurface");
memory->grow(nsurface, nmax, 3, "atom:rheo_nsurface");
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;
memory->grow(B, nmax, dim * dim, "rheo/surface:B");
memory->grow(gradC, nmax, dim * dim, "rheo/surface:gradC");
nmax_old = atom->nmax;
}
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;
B[i][a * dim + b] = 0.0;
gradC[i][a * dim + b] = 0.0;
}
n_surface[i][a] = 0.0;
nsurface[i][a] = 0.0;
}
divr[i] = 0.0;
surface[i] = 0;
// Remove surface settings
status[i] &= FixRHEO::surfacemask;
}
// loop over neighbors to calculate the average orientation of neighbors
@ -224,98 +182,105 @@ void ComputeRHEOSurface::pre_force(int /*vflag*/)
jlist = firstneigh[i];
jnum = numneigh[i];
itype = type[i];
fluidi = status[i] & FixRHEO::STATUS_FLUID;
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;
dx[0] = xtmp - x[j][0];
dx[1] = ytmp - x[j][1];
dx[2] = ztmp - x[j][2];
rsq = lensq(dx);
if (rsq < cutsq) {
jtype = type[j];
fluidj = status[j] & FixRHEO::STATUS_FLUID;
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;
if (fluidi && (!fluidj)) {
rhoj = compute_solids->correct_rho(j, i);
} else if ((!fluidi) && fluidj) {
rhoi = compute_solids->correct_rho(i, j);
} else if ((!fluidi) && (!fluidj)) {
rhoi = rho0;
rhoj = rho0;
}
Voli = mass[itype]/rhoi;
Volj = mass[jtype]/rhoj;
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);
wp = compute_kernel->calc_dw_quintic(i, j, dx[0], dx[1], dx[2], sqrt(rsq),dWij, 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;
for (a = 0; a < dim; a++){
divr[i] -= dWij[a] * dx[a] * Volj;
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;
for (a = 0; a < dim; a++){
divr[j] += dWji[a] * dx[a] * Voli;
gradC[j][a] += dWji[a] * Voli;
}
}
}
}
}
// reverse gradC and divr, forward divr
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);
comm_reverse = dim * dim + 1;
comm_forward = 1;
if (newton) comm->reverse_comm(this);
comm->forward_comm(this);
// calculate nsurface for local atoms
// Note, this isn't forwarded to ghosts
double maggC;
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
maggC = 0.0;
for (a = 0;a < dim; a++)
maggC += gradC[i][a] * gradC[i][a];
maggC = sqrt(maggC) + EPSILON;
maggC = 1.0 / maggC;
for (a = 0; a < dim; a++)
nsurface[i][a] = -gradC[i][a] * maggC;
}
}
int *coordination = compute_kernel->coordination;
// Find the free-surface
//0-bulk 1-surf vicinity 2-surface 3-splash
if (divr_flag) {
if (threshold_style == FixRHEO::DIVR) {
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;
status[i] |= FixRHEO::STATUS_BULK;
rsurface[i] = cut;
if (divr[i] < threshold_divr) {
status[i] |= FixRHEO::STATUS_SURFACE;
rsurface[i] = 0.0;
if (coordination[i] < threshold_z)
status[i] |= FixRHEO::STATUS_SPLASH;
}
}
}
} else {
for (i = 0; i < nall; i++) {
if (mask[i] & groupbit) {
surface[i] = 0;
rsurf[i] = cut; //Maximum range that can be seen
status[i] |= FixRHEO::STATUS_BULK;
rsurface[i] = cut;
if (coordination[i] < divR_limit) {
surface[i] = 2;
rsurf[i] = 0.0;
if (coordination[i] < coord_limit) surface[i] = 3;
status[i] |= FixRHEO::STATUS_SURFACE;
rsurface[i] = 0.0;
if (coordination[i] < threshold_z)
status[i] |= FixRHEO::STATUS_SPLASH;
}
}
}
}
//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];
@ -329,167 +294,37 @@ void ComputeRHEOSurface::pre_force(int /*vflag*/)
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;
dx[0] = xtmp - x[j][0];
dx[1] = ytmp - x[j][1];
dx[2] = ztmp - x[j][2];
rsq = lensq(dx);
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);
if ((status[i] & FixRHEO::STATUS_BULK) && (status[j] & FixRHEO::STATUS_SURFACE)) {
status[i] &= FixRHEO::surfacemask;
status[i] |= FixRHEO::STATUS_LAYER;
}
if (status[j] & FixRHEO::STATUS_SURFACE) rsurface[i] = MIN(rsurface[i], sqrt(rsq));
if (j < nlocal || newton) {
if ((status[j] & FixRHEO::STATUS_BULK) && (status[i] & FixRHEO::STATUS_SURFACE)) {
status[j] &= FixRHEO::surfacemask;
status[j] |= FixRHEO::STATUS_LAYER;
}
if (status[i] & FixRHEO::STATUS_SURFACE) rsurface[j] = MIN(rsurface[j], sqrt(rsq));
}
}
}
}
// forward/reverse status and rsurface
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
}
}
if (newton) comm->reverse_comm(this);
comm->forward_comm(this);
}
/* ---------------------------------------------------------------------- */
@ -498,8 +333,7 @@ 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];
int *status = atom->status;
m = 0;
last = first + n;
@ -508,14 +342,10 @@ int ComputeRHEOSurface::pack_reverse_comm(int n, int first, double *buf)
buf[m++] = divr[i];
for (a = 0; a < dim; a ++ )
for (b = 0; b < dim; b ++)
buf[m++] = gradC[i][a*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];
buf[m++] = (double) status[i];
buf[m++] = rsurface[i];
}
}
return m;
@ -527,8 +357,8 @@ 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];
int *status = atom->status;
int temp;
m = 0;
for (i = 0; i < n; i++) {
@ -537,16 +367,14 @@ void ComputeRHEOSurface::unpack_reverse_comm(int n, int *list, double *buf)
divr[j] += buf[m++];
for (a = 0; a < dim; a ++ )
for (b = 0; b < dim; b ++)
gradC[j][a*dim + b] += buf[m++];
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++];
temp = (int) buf[m++];
if ((status[j] & FixRHEO::STATUS_BULK) && (temp & FixRHEO::STATUS_LAYER))
status[j] = temp;
rsurface[j] = MIN(rsurface[j], buf[m++]);
}
}
}
@ -558,8 +386,7 @@ 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];
int *status = atom->status;
m = 0;
for (i = 0; i < n; i++) {
@ -567,8 +394,8 @@ int ComputeRHEOSurface::pack_forward_comm(int n, int *list, double *buf,
if (comm_stage == 0) {
buf[m++] = divr[j];
} else if (comm_stage == 1) {
buf[m++] = (double) surface[j];
buf[m++] = rsurf[j];
buf[m++] = (double) status[j];
buf[m++] = rsurface[j];
}
}
return m;
@ -579,8 +406,7 @@ int ComputeRHEOSurface::pack_forward_comm(int n, int *list, double *buf,
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];
int *status = atom->status;
m = 0;
last = first + n;
@ -588,8 +414,8 @@ void ComputeRHEOSurface::unpack_forward_comm(int n, int first, double *buf)
if (comm_stage == 0) {
divr[i] = buf[m++];
} else if (comm_stage == 1) {
surface[i] = (int) buf[m++];
rsurf[i] = buf[m++];
status[i] = (int) buf[m++];
rsurface[i] = buf[m++];
}
}
}

13
src/RHEO/compute_rheo_surface.h
View File

@ -36,16 +36,13 @@ class ComputeRHEOSurface : public Compute {
int pack_forward_comm(int, int *, double *, int, int *) override;
void unpack_forward_comm(int, int, double *) override;
double **gradC, **n_surface;
double **nsurface, **rsurface;
private:
double cut, cutsq, threshold;
int surface_style, nmax_old;
double **B, *divr;
int comm_stage;
int index_divr;
int index_rsurf;
double cut, cutsq, rho0, threshold_divr;
int surface_style, nmax_old, threshold_z;
double **B, **gradC, *divr;
int threshold_style, comm_stage;
double divR_limit;
int coord_limit;

2
src/RHEO/compute_rheo_vshift.cpp
View File

@ -36,7 +36,7 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeRHEOVShift::ComputeRHEOVShift(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), vshift(nullptr), fix_rheo(nullptr), fix_rheo(nullptr),
Compute(lmp, narg, arg), list(nullptr), vshift(nullptr), fix_rheo(nullptr),
compute_kernel(nullptr), compute_interface(nullptr)
{
if (narg != 3) error->all(FLERR,"Illegal compute RHEO/VShift command");

2
src/RHEO/fix_rheo.cpp
View File

@ -377,7 +377,7 @@ void FixRHEO::pre_force(int /*vflag*/)
if (shift_flag)
compute_vshift->compute_peratom();
// Remove extra shifting/no force options options
// Remove extra shifting/no force options
int *status = atom->status;
int nall = atom->nlocal + atom->nghost;
for (int i = 0; i < nall; i++) {

4
src/RHEO/fix_rheo_thermal.cpp
View File

@ -269,10 +269,10 @@ void FixRHEOThermal::post_integrate()
}
if (Ti > Tci) {
status[i] &= phasemask;
status[i] &= FixRHEO::phasemask;
status[i] |= FixRHEO::STATUS_FLUID;
} else if (!(status[i] & FixRHEO::STATUS_SOLID))
status[i] &= phasemask;
status[i] &= FixRHEO::phasemask;
status[i] |= FixRHEO::STATUS_FREEZING;
}
}

4
src/RHEO/pair_rheo.cpp
View File

@ -216,8 +216,8 @@ void PairRHEO::compute(int eflag, int vflag)
fmag = (chi[i] - 0.9) * (h * 0.5 - r) * rho0 * csq * h * rinv;
} else if ((!fluidi) && (!fluidj)) {
rhoi = 1.0;
rhoj = 1.0;
rhoi = rho0;
rhoj = rho0;
}
// Repel if close to inner solid particle