From a4d971df526e71f04a8de31b642ddd512100e89c Mon Sep 17 00:00:00 2001 From: jtclemm Date: Thu, 20 Apr 2023 14:45:35 -0600 Subject: [PATCH] Updating surface compute --- src/RHEO/compute_rheo_grad.cpp | 2 +- src/RHEO/compute_rheo_interface.cpp | 1 - src/RHEO/compute_rheo_kernel.cpp | 12 +- src/RHEO/compute_rheo_surface.cpp | 528 ++++++++++------------------ src/RHEO/compute_rheo_surface.h | 13 +- src/RHEO/compute_rheo_vshift.cpp | 2 +- src/RHEO/fix_rheo.cpp | 2 +- src/RHEO/fix_rheo_thermal.cpp | 4 +- src/RHEO/pair_rheo.cpp | 4 +- 9 files changed, 195 insertions(+), 373 deletions(-) diff --git a/src/RHEO/compute_rheo_grad.cpp b/src/RHEO/compute_rheo_grad.cpp index 31f60550c1..add3ed712d 100644 --- a/src/RHEO/compute_rheo_grad.cpp +++ b/src/RHEO/compute_rheo_grad.cpp @@ -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"); diff --git a/src/RHEO/compute_rheo_interface.cpp b/src/RHEO/compute_rheo_interface.cpp index 55cb81382b..059f4057f4 100644 --- a/src/RHEO/compute_rheo_interface.cpp +++ b/src/RHEO/compute_rheo_interface.cpp @@ -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); diff --git a/src/RHEO/compute_rheo_kernel.cpp b/src/RHEO/compute_rheo_kernel.cpp index ac9ea75fc4..2fd189e1e4 100644 --- a/src/RHEO/compute_rheo_kernel.cpp +++ b/src/RHEO/compute_rheo_kernel.cpp @@ -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; diff --git a/src/RHEO/compute_rheo_surface.cpp b/src/RHEO/compute_rheo_surface.cpp index dff5db4cfa..5b3aeabf26 100644 --- a/src/RHEO/compute_rheo_surface.cpp +++ b/src/RHEO/compute_rheo_surface.cpp @@ -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; areverse_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 rreverse_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;areverse_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++]; } } } diff --git a/src/RHEO/compute_rheo_surface.h b/src/RHEO/compute_rheo_surface.h index 480a8cd8ec..8d0b681c6f 100644 --- a/src/RHEO/compute_rheo_surface.h +++ b/src/RHEO/compute_rheo_surface.h @@ -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; diff --git a/src/RHEO/compute_rheo_vshift.cpp b/src/RHEO/compute_rheo_vshift.cpp index 6cbd6e96da..1b3edcffd8 100644 --- a/src/RHEO/compute_rheo_vshift.cpp +++ b/src/RHEO/compute_rheo_vshift.cpp @@ -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"); diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index 8fc48c0b3b..5358cc4cba 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -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++) { diff --git a/src/RHEO/fix_rheo_thermal.cpp b/src/RHEO/fix_rheo_thermal.cpp index 1912dc9f8c..fd08f39fd7 100644 --- a/src/RHEO/fix_rheo_thermal.cpp +++ b/src/RHEO/fix_rheo_thermal.cpp @@ -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; } } diff --git a/src/RHEO/pair_rheo.cpp b/src/RHEO/pair_rheo.cpp index 1eb2a43e1a..5974b9b756 100644 --- a/src/RHEO/pair_rheo.cpp +++ b/src/RHEO/pair_rheo.cpp @@ -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