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

Fixing compilation errors

Browse Source
This commit is contained in:
jtclemm
2023-04-24 19:46:27 -06:00
parent 35d1178cfa
commit 47b8cdc94f
21 changed files with 379 additions and 346 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/RHEO/atom_vec_rheo.cpp
View File

@ -120,7 +120,7 @@ void AtomVecRHEO::pack_property_atom(int index, double *buf, int nvalues, int gr
buf[n] = 0.0; buf[n] = 0.0;
n += nvalues; n += nvalues;
} }
} if else (index == 1) { } else if (index == 1) {
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) if (mask[i] & groupbit)
buf[n] = rho[i]; buf[n] = rho[i];

66
src/RHEO/compute_rheo_grad.cpp
View File

@ -21,11 +21,12 @@
#include "atom.h" #include "atom.h"
#include "comm.h" #include "comm.h"
#include "compute_rheo_kernel.h" #include "compute_rheo_kernel.h"
#include "compute_rheo_solids.h" #include "compute_rheo_interface.h"
#include "domain.h" #include "domain.h"
#include "error.h" #include "error.h"
#include "fix_rheo.h" #include "fix_rheo.h"
#include "force.h" #include "force.h"
#include "memory.h"
#include "neighbor.h" #include "neighbor.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "update.h" #include "update.h"
@ -33,6 +34,8 @@
#include <cmath> #include <cmath>
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace RHEO_NS;
enum{COMMGRAD, COMMFIELD}; enum{COMMGRAD, COMMFIELD};
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -112,10 +115,13 @@ void ComputeRHEOGrad::init()
compute_kernel = fix_rheo->compute_kernel; compute_kernel = fix_rheo->compute_kernel;
compute_interface = fix_rheo->compute_interface; compute_interface = fix_rheo->compute_interface;
int tmp1, tmp2;
index_visc = atom->find_custom("rheo_viscosity", tmp1, tmp2);
// Create coordination array if it doesn't already exist // Create coordination array if it doesn't already exist
// Create a custom atom property so it works with compute property/atom // 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 // Do not create grow callback as there's no reason to copy/exchange data
// Manually grow if nmax_old exceeded // Manually grow if nmax_store exceeded
int index; int index;
int dim = domain->dimension; int dim = domain->dimension;
@ -139,7 +145,7 @@ void ComputeRHEOGrad::init()
gradn = atom->darray[index]; gradn = atom->darray[index];
} }
nmax_old = 0; nmax_store = 0;
grow_arrays(atom->nmax); grow_arrays(atom->nmax);
} }
@ -158,7 +164,7 @@ void ComputeRHEOGrad::compute_peratom()
double xtmp, ytmp, ztmp, delx, dely, delz; double xtmp, ytmp, ztmp, delx, dely, delz;
double rsq, imass, jmass; double rsq, imass, jmass;
double rhoi, rhoj, Voli, Volj, drho, dT, deta; double rhoi, rhoj, Voli, Volj, drho, dT, deta;
double vij[3]; double vi[3], vj[3], vij[3];
double wp, *dWij, *dWji; double wp, *dWij, *dWji;
int inum, *ilist, *numneigh, **firstneigh; int inum, *ilist, *numneigh, **firstneigh;
@ -169,26 +175,22 @@ void ComputeRHEOGrad::compute_peratom()
double **v = atom->v; double **v = atom->v;
double *rho = atom->rho; double *rho = atom->rho;
double *temperature = atom->temperature; double *temperature = atom->temperature;
double *viscosity = atom->dvector[index_visc];
int *status = atom->status; int *status = atom->status;
int *type = atom->type; int *type = atom->type;
double *mass = atom->mass; double *mass = atom->mass;
int newton = force->newton; int newton = force->newton;
int dim = domain->dimension; int dim = domain->dimension;
int tmp1, tmp2;
int index_visc = atom->find_custom("rheo_viscosity", tmp1, tmp2);
if (index_visc == -1) error->all(FLERR, "Cannot find rheo viscosity");
double *viscosity = atom->dvector[index_visc];
inum = list->inum; inum = list->inum;
ilist = list->ilist; ilist = list->ilist;
numneigh = list->numneigh; numneigh = list->numneigh;
firstneigh = list->firstneigh; firstneigh = list->firstneigh;
// initialize arrays // initialize arrays
if (nmax > nmax_old) grow_arrays(nmax); if (atom->nmax > nmax_store) grow_arrays(atom->nmax);
for (i = 0; i < nmax; i++) { for (i = 0; i < nmax_store; i++) {
if (velocity_flag) { if (velocity_flag) {
for (k = 0; k < dim * dim; k++) for (k = 0; k < dim * dim; k++)
gradv[i][k] = 0.0; gradv[i][k] = 0.0;
@ -212,6 +214,9 @@ void ComputeRHEOGrad::compute_peratom()
xtmp = x[i][0]; xtmp = x[i][0];
ytmp = x[i][1]; ytmp = x[i][1];
ztmp = x[i][2]; ztmp = x[i][2];
vi[0] = v[i][0];
vi[1] = v[i][1];
vi[2] = v[i][2];
itype = type[i]; itype = type[i];
jlist = firstneigh[i]; jlist = firstneigh[i];
jnum = numneigh[i]; jnum = numneigh[i];
@ -230,14 +235,18 @@ void ComputeRHEOGrad::compute_peratom()
rhoi = rho[i]; rhoi = rho[i];
rhoj = rho[j]; rhoj = rho[j];
vj[0] = v[j][0];
vj[1] = v[j][1];
vj[2] = v[j][2];
// Add corrections for walls // Add corrections for walls
if ((status[i] & FixRHEO::STATUS_FLUID) && !(status[j] & FixRHEO::STATUS_FLUID)) { if ((status[i] & STATUS_FLUID) && !(status[j] & STATUS_FLUID)) {
compute_interface->correct_v(v[i], v[j], vi, i, j); compute_interface->correct_v(vi, vj, i, j);
rhoj = compute_interface->correct_rho(j,i); rhoj = compute_interface->correct_rho(j, i);
} else if (!(status[i] & FixRHEO::STATUS_FLUID) && (status[j] & FixRHEO::STATUS_FLUID)) { } else if (!(status[i] & STATUS_FLUID) && (status[j] & STATUS_FLUID)) {
compute_interface->correct_v(v[j], v[i], vj, j, i); compute_interface->correct_v(vj, vi, j, i);
rhoi = compute_interface->correct_rho(i,j); rhoi = compute_interface->correct_rho(i, j);
} else if (!(status[i] & FixRHEO::STATUS_FLUID) && !(status[j] & FixRHEO::STATUS_FLUID)) { } else if (!(status[i] & STATUS_FLUID) && !(status[j] & STATUS_FLUID)) {
rhoi = rho0; rhoi = rho0;
rhoj = rho0; rhoj = rho0;
} }
@ -324,7 +333,6 @@ int ComputeRHEOGrad::pack_forward_comm(int n, int *list, double *buf,
int i,j,k,m; int i,j,k,m;
double *rho = atom->rho; double *rho = atom->rho;
double *temperature = atom->temperature; double *temperature = atom->temperature;
double *eta = atom->viscosity;
double **v = atom->v; double **v = atom->v;
int dim = domain->dimension; int dim = domain->dimension;
@ -371,9 +379,9 @@ int ComputeRHEOGrad::pack_forward_comm(int n, int *list, double *buf,
void ComputeRHEOGrad::unpack_forward_comm(int n, int first, double *buf) void ComputeRHEOGrad::unpack_forward_comm(int n, int first, double *buf)
{ {
int i, k, m, last; int i, k, m, last;
double * rho = atom->rho; double *rho = atom->rho;
double * temperature = atom->temperature; double *temperature = atom->temperature;
double ** v = atom->v; double **v = atom->v;
int dim = domain->dimension; int dim = domain->dimension;
m = 0; m = 0;
@ -483,25 +491,27 @@ void ComputeRHEOGrad::grow_arrays(int nmax)
if (eta_flag) if (eta_flag)
memory->grow(gradn, nmax, dim, "atom:rheo_grad_eta"); memory->grow(gradn, nmax, dim, "atom:rheo_grad_eta");
nmax_old = nmax; nmax_store = nmax;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
double ComputeRHEOKernel::memory_usage() double ComputeRHEOGrad::memory_usage()
{ {
double bytes = 0.0; double bytes = 0.0;
int dim = domain->dimension;
if (velocity_flag) if (velocity_flag)
bytes = (size_t) nmax_old * dim * dim * sizeof(double); bytes = (size_t) nmax_store * dim * dim * sizeof(double);
if (rho_flag) if (rho_flag)
bytes = (size_t) nmax_old * dim * sizeof(double); bytes = (size_t) nmax_store * dim * sizeof(double);
if (temperature_flag) if (temperature_flag)
bytes = (size_t) nmax_old * dim * sizeof(double); bytes = (size_t) nmax_store * dim * sizeof(double);
if (eta_flag) if (eta_flag)
bytes = (size_t) nmax_old * dim * sizeof(double); bytes = (size_t) nmax_store * dim * sizeof(double);
return bytes; return bytes;
} }

5
src/RHEO/compute_rheo_grad.h
View File

@ -45,7 +45,8 @@ class ComputeRHEOGrad : public Compute {
class FixRHEO *fix_rheo; class FixRHEO *fix_rheo;
private: private:
int comm_stage, ncomm_grad, ncomm_field, nmax_old; int comm_stage, ncomm_grad, ncomm_field, nmax_store;
int index_visc;
double cut, cutsq, rho0; double cut, cutsq, rho0;
class NeighList *list; class NeighList *list;
@ -53,6 +54,8 @@ class ComputeRHEOGrad : public Compute {
class ComputeRHEOInterface *compute_interface; class ComputeRHEOInterface *compute_interface;
int velocity_flag, temperature_flag, rho_flag, eta_flag; int velocity_flag, temperature_flag, rho_flag, eta_flag;
void grow_arrays(int);
}; };
} // namespace LAMMPS_NS } // namespace LAMMPS_NS

91
src/RHEO/compute_rheo_interface.cpp
View File

@ -34,18 +34,19 @@
#include <cmath> #include <cmath>
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace RHEO_NS;
#define EPSILON 1e-1 static constexpr double EPSILON = 1e-1;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
ComputeRHEOInterface::ComputeRHEOInterface(LAMMPS *lmp, int narg, char **arg) : ComputeRHEOInterface::ComputeRHEOInterface(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), fix_rheo(nullptr), compute_kernel(nullptr), fx_m_norm(nullptr), Compute(lmp, narg, arg), fix_rheo(nullptr), compute_kernel(nullptr), fp_store(nullptr),
norm(nullptr), normwf(nullptr), chi(nullptr), f_pressure(nullptr), id_fix_pa(nullptr) norm(nullptr), normwf(nullptr), chi(nullptr), id_fix_pa(nullptr)
{ {
if (narg != 3) error->all(FLERR,"Illegal compute rheo/interface command"); if (narg != 3) error->all(FLERR,"Illegal compute rheo/interface command");
nmax = 0; nmax_store = 0;
comm_forward = 3; comm_forward = 3;
comm_reverse = 4; comm_reverse = 4;
@ -74,15 +75,15 @@ void ComputeRHEOInterface::init()
compute_kernel = fix_rheo->compute_kernel; compute_kernel = fix_rheo->compute_kernel;
rho0 = fix_rheo->rho0; rho0 = fix_rheo->rho0;
cut = fix_rheo->cut; cut = fix_rheo->cut;
cs = fix_rheo->cs; csq = fix_rheo->csq;
cs_inv = 1.0 / cs; csq_inv = 1.0 / csq;
cutsq = cut * cut; cutsq = cut * cut;
wall_max = sqrt(3.0) / 12.0 * cut; wall_max = sqrt(3.0) / 12.0 * cut;
// Create chi array if it doesn't already exist // Create chi array if it doesn't already exist
// Create a custom atom property so it works with compute property/atom // 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 // Do not create grow callback as there's no reason to copy/exchange data
// Manually grow if nmax_old exceeded // Manually grow if nmax_store exceeded
int tmp1, tmp2; int tmp1, tmp2;
int nmax = atom->nmax; int nmax = atom->nmax;
@ -93,7 +94,7 @@ void ComputeRHEOInterface::init()
memory->destroy(normwf); memory->destroy(normwf);
memory->create(norm, nmax, "rheo/interface:norm"); memory->create(norm, nmax, "rheo/interface:norm");
memory->create(normwf, nmax, "rheo/interface:normwf"); memory->create(normwf, nmax, "rheo/interface:normwf");
nmax_old = nmax; nmax_store = nmax;
} }
chi = atom->dvector[index]; chi = atom->dvector[index];
@ -104,13 +105,13 @@ void ComputeRHEOInterface::init()
index = atom->find_custom("fp_store", tmp1, tmp2); index = atom->find_custom("fp_store", tmp1, tmp2);
if (index == -1) { if (index == -1) {
id_fix_pa = utils::strdup(id + std::string("_fix_property_atom")); id_fix_pa = utils::strdup(id + std::string("_fix_property_atom"));
modify->add_fix(fmt::format("{} all property/atom d2_fp_store 3", id_fix_pa))); modify->add_fix(fmt::format("{} all property/atom d2_fp_store 3", id_fix_pa));
index = atom->find_custom("fp_store", tmp1, tmp2); index = atom->find_custom("fp_store", tmp1, tmp2);
} }
fp_store = atom->darray[index]; fp_store = atom->darray[index];
// need an occasional half neighbor list // need an occasional half neighbor list
neighbor->add_request(this, NeighConst::REQ_HALF); neighbor->add_request(this, NeighConst::REQ_DEFAULT);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -142,17 +143,17 @@ void ComputeRHEOInterface::compute_peratom()
numneigh = list->numneigh; numneigh = list->numneigh;
firstneigh = list->firstneigh; firstneigh = list->firstneigh;
if (atom->nmax > nmax_old) { if (atom->nmax > nmax_store) {
nmax_old = atom->nmax; nmax_store = atom->nmax;
memory->destroy(norm); memory->destroy(norm);
memory->destroy(normwf); memory->destroy(normwf);
memory->create(norm, nmax_old, "rheo/interface:norm"); memory->create(norm, nmax_store, "rheo/interface:norm");
memory->create(normwf, nmax_old, "rheo/interface:normwf"); memory->create(normwf, nmax_store, "rheo/interface:normwf");
memory->grow(chi, nmax_old, "rheo/interface:chi"); memory->grow(chi, nmax_store, "rheo/interface:chi");
} }
for (i = 0; i < nall; i++) { for (i = 0; i < nall; i++) {
if (!(status[i] & FixRHEO::STATUS_FLUID)) rho[i] = 0.0; if (!(status[i] & STATUS_FLUID)) rho[i] = 0.0;
normwf[i] = 0.0; normwf[i] = 0.0;
norm[i] = 0.0; norm[i] = 0.0;
chi[i] = 0.0; chi[i] = 0.0;
@ -164,7 +165,7 @@ void ComputeRHEOInterface::compute_peratom()
ytmp = x[i][1]; ytmp = x[i][1];
ztmp = x[i][2]; ztmp = x[i][2];
itype = type[i]; itype = type[i];
fluidi = status[i] & FixRHEO::STATUS_FLUID; fluidi = status[i] & STATUS_FLUID;
jlist = firstneigh[i]; jlist = firstneigh[i];
jnum = numneigh[i]; jnum = numneigh[i];
@ -179,12 +180,12 @@ void ComputeRHEOInterface::compute_peratom()
if (rsq < cutsq) { if (rsq < cutsq) {
jtype = type[j]; jtype = type[j];
fluidj = status[j] & FixRHEO::STATUS_FLUID; fluidj = status[j] & STATUS_FLUID;
w = compute_kernel->calc_w_quintic(i, j, delx, dely, delz, sqrt(rsq)); w = compute_kernel->calc_w_quintic(i, j, delx, dely, delz, sqrt(rsq));
status_match = 0; status_match = 0;
norm[i] += w; norm[i] += w;
if ((fluidi && fluidj) || ((!fluid) && (!fluidj))) if ((fluidi && fluidj) || ((!fluidi) && (!fluidj)))
status_match = 1; status_match = 1;
if (status_match) { if (status_match) {
@ -194,7 +195,7 @@ void ComputeRHEOInterface::compute_peratom()
dot = (-fp_store[0][j] + fp_store[0][i]) * delx; dot = (-fp_store[0][j] + fp_store[0][i]) * delx;
dot += (-fp_store[1][j] + fp_store[1][i]) * dely; dot += (-fp_store[1][j] + fp_store[1][i]) * dely;
dot += (-fp_store[2][j] + fp_store[2][i]) * delz; dot += (-fp_store[2][j] + fp_store[2][i]) * delz;
rho[i] += w * (cs * (rho[j] - rho0) - rho[j] * dot); rho[i] += w * (csq * (rho[j] - rho0) - rho[j] * dot);
normwf[i] += w; normwf[i] += w;
} }
} }
@ -208,7 +209,7 @@ void ComputeRHEOInterface::compute_peratom()
dot = (-fp_store[0][i] + fp_store[0][j]) * delx; dot = (-fp_store[0][i] + fp_store[0][j]) * delx;
dot += (-fp_store[1][i] + fp_store[1][j]) * dely; dot += (-fp_store[1][i] + fp_store[1][j]) * dely;
dot += (-fp_store[2][i] + fp_store[2][j]) * delz; dot += (-fp_store[2][i] + fp_store[2][j]) * delz;
rho[j] += w * (cs * (rho[i] - rho0) + rho[i] * dot); rho[j] += w * (csq * (rho[i] - rho0) + rho[i] * dot);
normwf[j] += w; normwf[j] += w;
} }
} }
@ -223,10 +224,10 @@ void ComputeRHEOInterface::compute_peratom()
if (norm[i] != 0.0) chi[i] /= norm[i]; if (norm[i] != 0.0) chi[i] /= norm[i];
// Recalculate rho for non-fluid particles // Recalculate rho for non-fluid particles
if (!(status[i] & FixRHEO::STATUS_FLUID)) { if (!(status[i] & STATUS_FLUID)) {
if (normwf[i] != 0.0) { if (normwf[i] != 0.0) {
// Stores rho for solid particles 1+Pw in Adami Adams 2012 // Stores rho for solid particles 1+Pw in Adami Adams 2012
rho[i] = MAX(EPSILON, rho0 + (rho[i] / normwf[i]) * cs_inv); rho[i] = MAX(EPSILON, rho0 + (rho[i] / normwf[i]) * csq_inv);
} else { } else {
rho[i] = rho0; rho[i] = rho0;
} }
@ -310,7 +311,7 @@ void ComputeRHEOInterface::unpack_reverse_comm(int n, int *list, double *buf)
j = list[i]; j = list[i];
norm[j] += buf[m++]; norm[j] += buf[m++];
chi[j] += buf[m++]; chi[j] += buf[m++];
if (!(status[j] & FixRHEO::STATUS_FLUID)){ if (!(status[j] & STATUS_FLUID)){
normwf[j] += buf[m++]; normwf[j] += buf[m++];
rho[j] += buf[m++]; rho[j] += buf[m++];
} else { } else {
@ -322,7 +323,7 @@ void ComputeRHEOInterface::unpack_reverse_comm(int n, int *list, double *buf)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void ComputeRHEOInterface::correct_v(double *vi, double *vj, double *vi_out, int i, int j) void ComputeRHEOInterface::correct_v(double *vi, double *vj, int i, int j)
{ {
double wall_prefactor, wall_denom, wall_numer; double wall_prefactor, wall_denom, wall_numer;
@ -333,9 +334,9 @@ void ComputeRHEOInterface::correct_v(double *vi, double *vj, double *vi_out, int
wall_prefactor = wall_numer / wall_denom; wall_prefactor = wall_numer / wall_denom;
vi_out[0] = (vi[0] - vj[0]) * wall_prefactor + vi[0]; vi[0] = (vi[0] - vj[0]) * wall_prefactor + vi[0];
vi_out[1] = (vi[1] - vj[1]) * wall_prefactor + vi[1]; vi[1] = (vi[1] - vj[1]) * wall_prefactor + vi[1];
vi_out[2] = (vi[2] - vj[2]) * wall_prefactor + vi[2]; vi[2] = (vi[2] - vj[2]) * wall_prefactor + vi[2];
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -352,28 +353,30 @@ double ComputeRHEOInterface::correct_rho(int i, int j)
void ComputeRHEOInterface::store_forces() void ComputeRHEOInterface::store_forces()
{ {
double minv; double minv;
double mass = atom->mass; int *type = atom->type;
double type = atom->type;
double **f = atom->f;
int *mask = atom->mask; int *mask = atom->mask;
double *mass = atom->mass;
double **f = atom->f;
// When this is called, fp_store stores the pressure force // When this is called, fp_store stores the pressure force
// After this method, fp_store instead stores non-pressure forces // After this method, fp_store instead stores non-pressure forces
// and is also normalized by the particles mass // and is also normalized by the particles mass
// If forces are overwritten by a fix, there are no pressure forces // If forces are overwritten by a fix, there are no pressure forces
// so just normalize // so just normalize
int ifix = modify->find_fix_by_style("setforce"); auto fixlist = modify->get_fix_by_style("setforce");
if (ifix != -1) { if (fixlist.size() == 0) {
for (int i = 0; i < atom->nlocal; i++) { for (const auto &fix : fixlist) {
minv = 1.0 / mass[type[i]]; for (int i = 0; i < atom->nlocal; i++) {
if (mask[i] & modify->fix[ifix]->groupbit) { minv = 1.0 / mass[type[i]];
fp_store[i][0] = f[i][0] * minv; if (mask[i] & fix->groupbit) {
fp_store[i][1] = f[i][1] * minv; fp_store[i][0] = f[i][0] * minv;
fp_store[i][2] = f[i][2] * minv; fp_store[i][1] = f[i][1] * minv;
} else { fp_store[i][2] = f[i][2] * minv;
fp_store[i][0] = (f[i][0] - fp_store[i][0]) * minv; } else {
fp_store[i][1] = (f[i][1] - fp_store[i][1]) * minv; fp_store[i][0] = (f[i][0] - fp_store[i][0]) * minv;
fp_store[i][2] = (f[i][2] - fp_store[i][2]) * minv; fp_store[i][1] = (f[i][1] - fp_store[i][1]) * minv;
fp_store[i][2] = (f[i][2] - fp_store[i][2]) * minv;
}
} }
} }
} else { } else {
@ -397,7 +400,7 @@ void ComputeRHEOInterface::store_forces()
double ComputeRHEOInterface::memory_usage() double ComputeRHEOInterface::memory_usage()
{ {
double bytes = 3 * nmax_old * sizeof(double); double bytes = 3 * nmax_store * sizeof(double);
return bytes; return bytes;
} }

10
src/RHEO/compute_rheo_interface.h
View File

@ -36,17 +36,17 @@ class ComputeRHEOInterface : public Compute {
int pack_reverse_comm(int, int, double *) override; int pack_reverse_comm(int, int, double *) override;
void unpack_reverse_comm(int, int *, double *) override; void unpack_reverse_comm(int, int *, double *) override;
double memory_usage() override; double memory_usage() override;
void correct_v(double *, double *, double *, int, int); void correct_v(double *, double *, int, int);
double correct_rho(int, int); double correct_rho(int, int);
void store_forces(); void store_forces();
double *chi, **f_pressure; double *chi, **fp_store;
class FixRHEO *fix_rheo; class FixRHEO *fix_rheo;
private: private:
int nmax_old, comm_stage; int nmax_store, comm_stage;
double rho0, cut, cutsq, cs, cs_inv, wall_max; double rho0, cut, cutsq, csq, csq_inv, wall_max;
double *norm, *normwf, **fom_store; double *norm, *normwf;
char *id_fix_pa; char *id_fix_pa;

112
src/RHEO/compute_rheo_kernel.cpp
View File

@ -42,10 +42,10 @@
#include <gsl/gsl_cblas.h> #include <gsl/gsl_cblas.h>
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace RHEO_NS;
using namespace MathExtra; using namespace MathExtra;
enum {QUINTIC, CRK0, CRK1, CRK2}; static constexpr int DELTA = 2000;
#define DELTA 2000
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -55,16 +55,16 @@ ComputeRHEOKernel::ComputeRHEOKernel(LAMMPS *lmp, int narg, char **arg) :
{ {
if (narg != 4) error->all(FLERR,"Illegal compute rheo/kernel command"); if (narg != 4) error->all(FLERR,"Illegal compute rheo/kernel command");
kernel_style = (SubModelType) utils::inumeric(FLERR,arg[3],false,lmp); kernel_style = utils::inumeric(FLERR,arg[3],false,lmp);
if (kernel_style == FixRHEO::QUINTIC) { if (kernel_style == QUINTIC) {
correction_order = -1; correction_order = -1;
} else if (kernel_style == FixRHEO::CRK0) { } else if (kernel_style == CRK0) {
correction_order = 0; correction_order = 0;
} else if (kernel_style == FixRHEO::CRK1) { } else if (kernel_style == CRK1) {
correction_order = 1; correction_order = 1;
} else if (kernel_style == FixRHEO::CRK2) { } else if (kernel_style == CRK2) {
correction_order = 2; correction_order = 2;
} }
@ -74,11 +74,11 @@ ComputeRHEOKernel::ComputeRHEOKernel(LAMMPS *lmp, int narg, char **arg) :
comm_forward = 1; comm_forward = 1;
ncor = 0; ncor = 0;
Mdim = 0; Mdim = 0;
if (kernel_type == CRK1) { if (kernel_style == CRK1) {
Mdim = 1 + dim; Mdim = 1 + dim;
ncor = 1 + dim; ncor = 1 + dim;
comm_forward = ncor * Mdim; comm_forward = ncor * Mdim;
} else if (kernel_type == CRK2) { } else if (kernel_style == CRK2) {
//Polynomial basis size (up to quadratic order) //Polynomial basis size (up to quadratic order)
Mdim = 1 + dim + dim * (dim + 1) / 2; Mdim = 1 + dim + dim * (dim + 1) / 2;
//Number of sets of correction coefficients (1 x y xx yy) + z zz (3D) //Number of sets of correction coefficients (1 x y xx yy) + z zz (3D)
@ -126,35 +126,35 @@ void ComputeRHEOKernel::init()
if (dim == 3) { if (dim == 3) {
pre_w = 0.002652582384864922 * 27.0 * ihsq * ih; pre_w = 0.002652582384864922 * 27.0 * hsqinv * hinv;
pre_wp = pre_w * 3.0 * ih; pre_wp = pre_w * 3.0 * hinv;
} else { } else {
pre_w = 0.004661441847879780 * 9 * ihsq; pre_w = 0.004661441847879780 * 9 * hsqinv;
pre_wp = pre_w * 3.0 * ih; pre_wp = pre_w * 3.0 * hinv;
} }
// Create coordination array if it doesn't already exist // Create coordination array if it doesn't already exist
// Create a custom atom property so it works with compute property/atom // 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 // Do not create grow callback as there's no reason to copy/exchange data
// Manually grow if nmax_old exceeded // Manually grow if nmax_store exceeded
int tmp1, tmp2; int tmp1, tmp2;
int nmax = atom->nmax; int nmax = atom->nmax;
int index = atom->find_custom("rheo_coordination", tmp1, tmp2); int index = atom->find_custom("rheo_coordination", tmp1, tmp2);
if (index == -1) { if (index == -1) {
index = atom->add_custom("rheo_coordination", 0, 0); index = atom->add_custom("rheo_coordination", 0, 0);
nmax_old = nmax; nmax_store = nmax;
} }
coordination = atom->ivector[index]; coordination = atom->ivector[index];
// Create local arrays for kernel arrays, I can't foresee a reason to print // Create local arrays for kernel arrays, I can't foresee a reason to print
if (kernel_type == CRK0) { if (kernel_style == CRK0) {
memory->create(C0, nmax_old, "rheo/kernel:C0"); memory->create(C0, nmax_store, "rheo/kernel:C0");
} else if (kernel_type == CRK1) { } else if (kernel_style == CRK1) {
memory->create(C, nmax_old, ncor, Mdim, "rheo/kernel:C"); memory->create(C, nmax_store, ncor, Mdim, "rheo/kernel:C");
} else if (kernel_type == CRK2) { } else if (kernel_style == CRK2) {
memory->create(C, nmax_old, ncor, Mdim, "rheo/kernel:C"); memory->create(C, nmax_store, ncor, Mdim, "rheo/kernel:C");
} }
} }
@ -192,10 +192,10 @@ double ComputeRHEOKernel::calc_w(int i, int j, double delx, double dely, double
int corrections_j = check_corrections(j); int corrections_j = check_corrections(j);
int corrections = corrections_i & corrections_j; int corrections = corrections_i & corrections_j;
if (kernel_type == QUINTIC || !corrections) w = calc_w_quintic(i,j,delx,dely,delz,r); if (kernel_style == QUINTIC || !corrections) w = calc_w_quintic(i,j,delx,dely,delz,r);
else if (kernel_type == CRK0) w = calc_w_crk0(i,j,delx,dely,delz,r); else if (kernel_style == CRK0) w = calc_w_crk0(i,j,delx,dely,delz,r);
else if (kernel_type == CRK1) w = calc_w_crk1(i,j,delx,dely,delz,r); else if (kernel_style == CRK1) w = calc_w_crk1(i,j,delx,dely,delz,r);
else if (kernel_type == CRK2) w = calc_w_crk2(i,j,delx,dely,delz,r); else if (kernel_style == CRK2) w = calc_w_crk2(i,j,delx,dely,delz,r);
return w; return w;
} }
@ -211,11 +211,11 @@ double ComputeRHEOKernel::calc_dw(int i, int j, double delx, double dely, double
// Calc wp and default dW's, a bit inefficient but can redo later // Calc wp and default dW's, a bit inefficient but can redo later
wp = calc_dw_quintic(i,j,delx,dely,delz,r,dWij,dWji); wp = calc_dw_quintic(i,j,delx,dely,delz,r,dWij,dWji);
if(kernel_type == CRK1) { if(kernel_style == CRK1) {
//check if kernel correction calculated successfully. If not, revert to quintic //check if kernel correction calculated successfully. If not, revert to quintic
if (corrections_i) calc_dw_crk1(i,j,delx,dely,delz,r,dWij); if (corrections_i) calc_dw_crk1(i,j,delx,dely,delz,r,dWij);
if (corrections_j) calc_dw_crk1(j,i,-delx,-dely,-delz,r,dWji); if (corrections_j) calc_dw_crk1(j,i,-delx,-dely,-delz,r,dWji);
} else if(kernel_type == CRK2) { } else if(kernel_style == CRK2) {
if (corrections_i) calc_dw_crk2(i,j,delx,dely,delz,r,dWij); if (corrections_i) calc_dw_crk2(i,j,delx,dely,delz,r,dWij);
if (corrections_j) calc_dw_crk2(j,i,-delx,-dely,-delz,r,dWji); if (corrections_j) calc_dw_crk2(j,i,-delx,-dely,-delz,r,dWji);
} }
@ -228,7 +228,7 @@ double ComputeRHEOKernel::calc_dw(int i, int j, double delx, double dely, double
double ComputeRHEOKernel::calc_w_quintic(int i, int j, double delx, double dely, double delz, double r) double ComputeRHEOKernel::calc_w_quintic(int i, int j, double delx, double dely, double delz, double r)
{ {
double w, tmp1, tmp2, tmp3, tmp1sq, tmp2sq, tmp3sq, s; double w, tmp1, tmp2, tmp3, tmp1sq, tmp2sq, tmp3sq, s;
s = r * 3.0 * ih; s = r * 3.0 * hinv;
if (s > 3.0) { if (s > 3.0) {
w = 0.0; w = 0.0;
@ -266,7 +266,7 @@ double ComputeRHEOKernel::calc_dw_quintic(int i, int j, double delx, double dely
double *mass = atom->mass; double *mass = atom->mass;
int *type = atom->type; int *type = atom->type;
s = r * 3.0 * ih; s = r * 3.0 * hinv;
if (s > 3.0) { if (s > 3.0) {
wp = 0.0; wp = 0.0;
@ -496,9 +496,9 @@ void ComputeRHEOKernel::compute_peratom()
gsl_error_flag = 0; gsl_error_flag = 0;
gsl_error_tags.clear(); gsl_error_tags.clear();
if (kernel_type == FixRHEO::QUINTIC) return; if (kernel_style == QUINTIC) return;
int i, j, ii, jj, inum, jnum, g, a, b, gsl_error; int i, j, ii, jj, inum, jnum, itype, g, a, b, gsl_error;
double xtmp, ytmp, ztmp, r, rsq, w, vj; double xtmp, ytmp, ztmp, r, rsq, w, vj;
double dx[3]; double dx[3];
gsl_matrix_view gM; gsl_matrix_view gM;
@ -520,9 +520,9 @@ void ComputeRHEOKernel::compute_peratom()
firstneigh = list->firstneigh; firstneigh = list->firstneigh;
// Grow arrays if necessary // Grow arrays if necessary
if (nmax_old < atom->nmax) grow_arrays(atom->nmax); if (nmax_store < atom->nmax) grow_arrays(atom->nmax);
if (kernel_type == FixRHEO::CRK0) { if (kernel_style == CRK0) {
double M; double M;
for (ii = 0; ii < inum; ii++) { for (ii = 0; ii < inum; ii++) {
@ -544,12 +544,12 @@ void ComputeRHEOKernel::compute_peratom()
dx[0] = xtmp - x[j][0]; dx[0] = xtmp - x[j][0];
dx[1] = ytmp - x[j][1]; dx[1] = ytmp - x[j][1];
dx[2] = ztmp - x[j][2]; dx[2] = ztmp - x[j][2];
rsq = lensq(dx); rsq = lensq3(dx);
if (rsq < hsq) { if (rsq < hsq) {
r = sqrt(rsq); r = sqrt(rsq);
w = calc_w_quintic(i,j,dx[0],dx[1],dx[2],r); w = calc_w_quintic(i,j,dx[0],dx[1],dx[2],r);
if (!(status[j] & FixRHEO::STATUS_FLUID) && solid_flag) { if (!(status[j] & STATUS_FLUID) && solid_flag) {
vj = mass[type[j]] / compute_interface->correct_rho(j,i); vj = mass[type[j]] / compute_interface->correct_rho(j,i);
} else vj = mass[type[j]] / rho[j]; } else vj = mass[type[j]] / rho[j];
@ -590,22 +590,24 @@ void ComputeRHEOKernel::compute_peratom()
dx[1] = ytmp - x[j][1]; dx[1] = ytmp - x[j][1];
dx[2] = ztmp - x[j][2]; dx[2] = ztmp - x[j][2];
rsq = lensq(dx); rsq = lensq3(dx);
if (rsq < cutsq) { if (rsq < hsq) {
r = sqrt(rsq); r = sqrt(rsq);
w = calc_w_quintic(i,j,dx[0],dx[1],dx[2],r); w = calc_w_quintic(i,j,dx[0],dx[1],dx[2],r);
if (status[j] > FixRHEO::FLUID_MAX && solid_flag) if (solid_flag)
vj = mass[type[j]]/compute_interface->correct_rho(j,i); if (!(status[j] & STATUS_FLUID))
else vj = mass[type[j]]/rho[j]; vj = mass[type[j]]/compute_interface->correct_rho(j,i);
else
vj = mass[type[j]]/rho[j];
//Populate the H-vector of polynomials (2D) //Populate the H-vector of polynomials (2D)
if (dim == 2) { if (dim == 2) {
H[0] = 1.0; H[0] = 1.0;
H[1] = dx[0] * hinv; H[1] = dx[0] * hinv;
H[2] = dx[1] * hinv; H[2] = dx[1] * hinv;
if (kernel_type == FixRHEO::CRK2) { if (kernel_style == CRK2) {
H[3] = 0.5 * dx[0] * dx[0] * hsqinv; H[3] = 0.5 * dx[0] * dx[0] * hsqinv;
H[4] = 0.5 * dx[1] * dx[1] * hsqinv; H[4] = 0.5 * dx[1] * dx[1] * hsqinv;
H[5] = dx[0] * dx[1] * hsqinv; H[5] = dx[0] * dx[1] * hsqinv;
@ -615,7 +617,7 @@ void ComputeRHEOKernel::compute_peratom()
H[1] = dx[0] * hinv; H[1] = dx[0] * hinv;
H[2] = dx[1] * hinv; H[2] = dx[1] * hinv;
H[3] = dx[2] * hinv; H[3] = dx[2] * hinv;
if (kernel_type == FixRHEO::CRK2) { if (kernel_style == CRK2) {
H[4] = 0.5 * dx[0] * dx[0] * hsqinv; H[4] = 0.5 * dx[0] * dx[0] * hsqinv;
H[5] = 0.5 * dx[1] * dx[1] * hsqinv; H[5] = 0.5 * dx[1] * dx[1] * hsqinv;
H[6] = 0.5 * dx[2] * dx[2] * hsqinv; H[6] = 0.5 * dx[2] * dx[2] * hsqinv;
@ -642,7 +644,7 @@ void ComputeRHEOKernel::compute_peratom()
} }
// Skip if undercoordinated // Skip if undercoordinated
if (coordination[i] < zmin) continue if (coordination[i] < zmin) continue;
// Use gsl to get Minv, use Cholesky decomposition since the // Use gsl to get Minv, use Cholesky decomposition since the
// polynomials are independent, M is symmetrix & positive-definite // polynomials are independent, M is symmetrix & positive-definite
@ -656,7 +658,7 @@ void ComputeRHEOKernel::compute_peratom()
//check if not positive-definite //check if not positive-definite
if (gsl_error != GSL_EDOM) if (gsl_error != GSL_EDOM)
error->warn(FLERR, "Failed decomposition in rheo/kernel, gsl_error = {}", gsl_error); error->warning(FLERR, "Failed decomposition in rheo/kernel, gsl_error = {}", gsl_error);
continue; continue;
} }
@ -689,7 +691,7 @@ void ComputeRHEOKernel::compute_peratom()
// columns 1-2 (2D) or 1-3 (3D) // columns 1-2 (2D) or 1-3 (3D)
//Second derivatives //Second derivatives
if (kernel_type == FixRHEO::CRK2) if (kernel_style == CRK2)
C[i][1 + dim + b][a] = M[a * Mdim + b + 1 + dim] * hsqinv; C[i][1 + dim + b][a] = M[a * Mdim + b + 1 + dim] * hsqinv;
// columns 3-4 (2D) or 4-6 (3D) // columns 3-4 (2D) or 4-6 (3D)
} }
@ -721,7 +723,7 @@ void ComputeRHEOKernel::compute_coordination()
firstneigh = list->firstneigh; firstneigh = list->firstneigh;
// Grow arrays if necessary // Grow arrays if necessary
if (nmax_old < atom->nmax) grow_arrays(atom->nmax); if (nmax_store < atom->nmax) grow_arrays(atom->nmax);
for (ii = 0; ii < inum; ii++) { for (ii = 0; ii < inum; ii++) {
i = ilist[ii]; i = ilist[ii];
@ -740,7 +742,7 @@ void ComputeRHEOKernel::compute_coordination()
dx[0] = xtmp - x[j][0]; dx[0] = xtmp - x[j][0];
dx[1] = ytmp - x[j][1]; dx[1] = ytmp - x[j][1];
dx[2] = ztmp - x[j][2]; dx[2] = ztmp - x[j][2];
rsq = lensq(dx); rsq = lensq3(dx);
if (rsq < hsq) if (rsq < hsq)
coordination[i] += 1; coordination[i] += 1;
@ -759,13 +761,13 @@ void ComputeRHEOKernel::grow_arrays(int nmax)
{ {
memory->grow(coordination, nmax, "atom:rheo_coordination"); memory->grow(coordination, nmax, "atom:rheo_coordination");
if (kernel_type == FixRHEO::CRK0) { if (kernel_style == CRK0) {
memory->grow(C0, nmax, "rheo/kernel:C0"); memory->grow(C0, nmax, "rheo/kernel:C0");
} else if (correction_order > 0) { } else if (correction_order > 0) {
memory->grow(C, nmax, ncor, Mdim, "rheo/kernel:C"); memory->grow(C, nmax, ncor, Mdim, "rheo/kernel:C");
} }
nmax_old = nmax; nmax_store = nmax;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -781,7 +783,7 @@ int ComputeRHEOKernel::pack_forward_comm(int n, int *list, double *buf,
if (comm_stage == 0) { if (comm_stage == 0) {
buf[m++] = coordination[j]; buf[m++] = coordination[j];
} else { } else {
if (kernel_type == FixRHEO::CRK0) { if (kernel_style == CRK0) {
buf[m++] = C0[j]; buf[m++] = C0[j];
} else { } else {
for (a = 0; a < ncor; a++) for (a = 0; a < ncor; a++)
@ -805,7 +807,7 @@ void ComputeRHEOKernel::unpack_forward_comm(int n, int first, double *buf)
if (comm_stage == 0) { if (comm_stage == 0) {
coordination[i] = buf[m++]; coordination[i] = buf[m++];
} else { } else {
if (kernel_type == FixRHEO::CRK0) { if (kernel_style == CRK0) {
C0[i] = buf[m++]; C0[i] = buf[m++];
} else { } else {
for (a = 0; a < ncor; a++) for (a = 0; a < ncor; a++)
@ -821,12 +823,12 @@ void ComputeRHEOKernel::unpack_forward_comm(int n, int first, double *buf)
double ComputeRHEOKernel::memory_usage() double ComputeRHEOKernel::memory_usage()
{ {
double bytes = 0.0; double bytes = 0.0;
bytes = (size_t) nmax_old * sizeof(int); bytes = (size_t) nmax_store * sizeof(int);
if (kernel_type == FixRHEO::CRK0) { if (kernel_style == CRK0) {
bytes += (size_t) nmax_old * sizeof(double); bytes += (size_t) nmax_store * sizeof(double);
} else if (correction_order > 0) { } else if (correction_order > 0) {
bytes += (size_t) nmax_old * ncor * Mdim * sizeof(double); bytes += (size_t) nmax_store * ncor * Mdim * sizeof(double);
} }
return bytes; return bytes;
} }

2
src/RHEO/compute_rheo_kernel.h
View File

@ -54,7 +54,7 @@ class ComputeRHEOKernel : public Compute {
std::unordered_set<tagint> gsl_error_tags; std::unordered_set<tagint> gsl_error_tags;
int kernel_style, zmin, dim, Mdim, ncor; int kernel_style, zmin, dim, Mdim, ncor;
int nmax_old; int nmax_store;
double h, hsq, hinv, hsqinv, pre_w, pre_wp; double h, hsq, hinv, hsqinv, pre_w, pre_wp;
double ***C; double ***C;
double *C0; double *C0;

13
src/RHEO/compute_rheo_rho_sum.cpp
View File

@ -49,7 +49,7 @@ void ComputeRHEORhoSum::init()
cutsq = cut * cut; cutsq = cut * cut;
// need an occasional half neighbor list // need an occasional half neighbor list
neighbor->add_request(this, NeighConst::REQ_HALF); neighbor->add_request(this, NeighConst::REQ_DEFAULT);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -130,8 +130,9 @@ void ComputeRHEORhoSum::compute_peratom()
int ComputeRHEORhoSum::pack_forward_comm(int n, int *list, double *buf, int ComputeRHEORhoSum::pack_forward_comm(int n, int *list, double *buf,
int /*pbc_flag*/, int * /*pbc*/) int /*pbc_flag*/, int * /*pbc*/)
{ {
int i,j,k,m; int i, j, k, m;
double * rho = atom->rho; double *rho = atom->rho;
int *coordination = compute_kernel->coordination;
m = 0; m = 0;
for (i = 0; i < n; i++) { for (i = 0; i < n; i++) {
@ -145,7 +146,7 @@ int ComputeRHEORhoSum::pack_forward_comm(int n, int *list, double *buf,
void ComputeRHEORhoSum::unpack_forward_comm(int n, int first, double *buf) void ComputeRHEORhoSum::unpack_forward_comm(int n, int first, double *buf)
{ {
int i, k, m, last; int i, k, m, last;
double * rho = atom->rho; double *rho = atom->rho;
m = 0; m = 0;
last = first + n; last = first + n;
@ -158,7 +159,7 @@ void ComputeRHEORhoSum::unpack_forward_comm(int n, int first, double *buf)
int ComputeRHEORhoSum::pack_reverse_comm(int n, int first, double *buf) int ComputeRHEORhoSum::pack_reverse_comm(int n, int first, double *buf)
{ {
int i,k,m,last; int i, k, m, last;
double *rho = atom->rho; double *rho = atom->rho;
m = 0; m = 0;
@ -173,7 +174,7 @@ int ComputeRHEORhoSum::pack_reverse_comm(int n, int first, double *buf)
void ComputeRHEORhoSum::unpack_reverse_comm(int n, int *list, double *buf) void ComputeRHEORhoSum::unpack_reverse_comm(int n, int *list, double *buf)
{ {
int i,k,j,m; int i, k, j, m;
double *rho = atom->rho; double *rho = atom->rho;
m = 0; m = 0;

83
src/RHEO/compute_rheo_surface.cpp
View File

@ -20,8 +20,8 @@
#include "atom.h" #include "atom.h"
#include "comm.h" #include "comm.h"
#include "compute_rheo_interface.h"
#include "compute_rheo_kernel.h" #include "compute_rheo_kernel.h"
#include "compute_rheo_solids.h"
#include "domain.h" #include "domain.h"
#include "error.h" #include "error.h"
#include "fix_rheo.h" #include "fix_rheo.h"
@ -33,18 +33,19 @@
#include "neigh_request.h" #include "neigh_request.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace RHEO_NS;
using namespace FixConst; using namespace FixConst;
using namespace MathExtra; using namespace MathExtra;
#define EPSILON 1e-10; static constexpr double EPSILON = 1e-10;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
ComputeRHEOSurface::ComputeRHEOSurface(LAMMPS *lmp, int narg, char **arg) : ComputeRHEOSurface::ComputeRHEOSurface(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), fix_rheo(nullptr), list(nullptr), compute_kernel(nullptr), compute_solids(nullptr), Compute(lmp, narg, arg), fix_rheo(nullptr), list(nullptr), compute_kernel(nullptr), compute_interface(nullptr),
B(nullptr), gradC(nullptr), nsurface(nullptr), divr(nullptr), rsurface(nullptr) B(nullptr), gradC(nullptr), nsurface(nullptr), divr(nullptr), rsurface(nullptr)
{ {
if (narg != 3) error->all(FLERR,"Illegal fix RHEO/SURFACE command"); if (narg != 3) error->all(FLERR,"Illegal compute RHEO/SURFACE command");
int dim = domain->dimension; int dim = domain->dimension;
comm_forward = 2; comm_forward = 2;
@ -75,19 +76,19 @@ ComputeRHEOSurface::~ComputeRHEOSurface()
void ComputeRHEOSurface::init() void ComputeRHEOSurface::init()
{ {
compute_kernel = fix_rheo->compute_kernel; compute_kernel = fix_rheo->compute_kernel;
compute_solids = fix_rheo->compute_solids; compute_interface = fix_rheo->compute_interface;
cut = fix_rheo->cut; cut = fix_rheo->cut;
rho0 = fix_rheo->rho0; rho0 = fix_rheo->rho0;
threshold_style = fix_rheo->surface_style; threshold_style = fix_rheo->surface_style;
threshold_divr = fix_rheo->divrsurface; threshold_divr = fix_rheo->divr_surface;
threshold_z = fix_rheo->zminsurface; threshold_z = fix_rheo->zmin_surface;
cutsq = cut * cut; cutsq = cut * cut;
// Create rsurface, divr, nsurface arrays if they don't already exist // Create rsurface, divr, nsurface arrays if they don't already exist
// Create a custom atom property so it works with compute property/atom // 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 // Do not create grow callback as there's no reason to copy/exchange data
// Manually grow if nmax_old exceeded // Manually grow if nmax_store exceeded
// For B and gradC, create a local array since they are unlikely to be printed // For B and gradC, create a local array since they are unlikely to be printed
int tmp1, tmp2; int tmp1, tmp2;
@ -103,12 +104,13 @@ void ComputeRHEOSurface::init()
if (index == -1) index = atom->add_custom("rheo_nsurface", 1, 3); if (index == -1) index = atom->add_custom("rheo_nsurface", 1, 3);
nsurface = atom->darray[index]; nsurface = atom->darray[index];
nmax_old = atom->nmax; nmax_store = atom->nmax;
memory->create(B, nmax_old, dim * dim, "rheo/surface:B"); int dim = domain->dimension;
memory->create(gradC, nmax_old, dim * dim, "rheo/surface:gradC"); memory->create(B, nmax_store, dim * dim, "rheo/surface:B");
memory->create(gradC, nmax_store, dim * dim, "rheo/surface:gradC");
// need an occasional half neighbor list // need an occasional half neighbor list
neighbor->add_request(this, NeighConst::REQ_HALF); neighbor->add_request(this, NeighConst::REQ_DEFAULT);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -123,9 +125,8 @@ void ComputeRHEOSurface::init_list(int /*id*/, NeighList *ptr)
void ComputeRHEOSurface::compute_peratom() void ComputeRHEOSurface::compute_peratom()
{ {
int i, j, ii, jj, inum, jnum, a, b, itype, jtype, fluidi, fluidj; int i, j, ii, jj, inum, jnum, a, b, itype, jtype, fluidi, fluidj;
double xtmp, ytmp, ztmp, rsq, Voli, Volj, rhoi, rhoj; double xtmp, ytmp, ztmp, rsq, Voli, Volj, rhoi, rhoj, wp;
double *dWij, *dWji; double *dWij, *dWji, dx[3];
double dx[3];
int *ilist, *jlist, *numneigh, **firstneigh; int *ilist, *jlist, *numneigh, **firstneigh;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
@ -146,7 +147,7 @@ void ComputeRHEOSurface::compute_peratom()
firstneigh = list->firstneigh; firstneigh = list->firstneigh;
int nmax = atom->nmax; int nmax = atom->nmax;
if (nmax_old <= nmax) { if (nmax_store <= nmax) {
memory->grow(divr, nmax, "atom:rheo_divr"); memory->grow(divr, nmax, "atom:rheo_divr");
memory->grow(rsurface, nmax, "atom:rheo_rsurface"); memory->grow(rsurface, nmax, "atom:rheo_rsurface");
memory->grow(nsurface, nmax, 3, "atom:rheo_nsurface"); memory->grow(nsurface, nmax, 3, "atom:rheo_nsurface");
@ -154,7 +155,7 @@ void ComputeRHEOSurface::compute_peratom()
memory->grow(B, nmax, dim * dim, "rheo/surface:B"); memory->grow(B, nmax, dim * dim, "rheo/surface:B");
memory->grow(gradC, nmax, dim * dim, "rheo/surface:gradC"); memory->grow(gradC, nmax, dim * dim, "rheo/surface:gradC");
nmax_old = atom->nmax; nmax_store = atom->nmax;
} }
int nall = nlocal + atom->nghost; int nall = nlocal + atom->nghost;
@ -169,7 +170,7 @@ void ComputeRHEOSurface::compute_peratom()
divr[i] = 0.0; divr[i] = 0.0;
// Remove surface settings // Remove surface settings
status[i] &= FixRHEO::surfacemask; status[i] &= SURFACEMASK;
} }
// loop over neighbors to calculate the average orientation of neighbors // loop over neighbors to calculate the average orientation of neighbors
@ -182,7 +183,7 @@ void ComputeRHEOSurface::compute_peratom()
jlist = firstneigh[i]; jlist = firstneigh[i];
jnum = numneigh[i]; jnum = numneigh[i];
itype = type[i]; itype = type[i];
fluidi = status[i] & FixRHEO::STATUS_FLUID; fluidi = status[i] & STATUS_FLUID;
for (jj = 0; jj < jnum; jj++) { for (jj = 0; jj < jnum; jj++) {
j = jlist[jj]; j = jlist[jj];
@ -192,19 +193,19 @@ void ComputeRHEOSurface::compute_peratom()
dx[1] = ytmp - x[j][1]; dx[1] = ytmp - x[j][1];
dx[2] = ztmp - x[j][2]; dx[2] = ztmp - x[j][2];
rsq = lensq(dx); rsq = lensq3(dx);
if (rsq < cutsq) { if (rsq < cutsq) {
jtype = type[j]; jtype = type[j];
fluidj = status[j] & FixRHEO::STATUS_FLUID; fluidj = status[j] & STATUS_FLUID;
rhoi = rho[i]; rhoi = rho[i];
rhoj = rho[j]; rhoj = rho[j];
// Add corrections for walls // Add corrections for walls
if (fluidi && (!fluidj)) { if (fluidi && (!fluidj)) {
rhoj = compute_solids->correct_rho(j, i); rhoj = compute_interface->correct_rho(j, i);
} else if ((!fluidi) && fluidj) { } else if ((!fluidi) && fluidj) {
rhoi = compute_solids->correct_rho(i, j); rhoi = compute_interface->correct_rho(i, j);
} else if ((!fluidi) && (!fluidj)) { } else if ((!fluidi) && (!fluidj)) {
rhoi = rho0; rhoi = rho0;
rhoj = rho0; rhoj = rho0;
@ -253,29 +254,29 @@ void ComputeRHEOSurface::compute_peratom()
} }
// Find the free-surface // Find the free-surface
if (threshold_style == FixRHEO::DIVR) { if (threshold_style == DIVR) {
for (i = 0; i < nall; i++) { for (i = 0; i < nall; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
status[i] |= FixRHEO::STATUS_BULK; status[i] |= STATUS_BULK;
rsurface[i] = cut; rsurface[i] = cut;
if (divr[i] < threshold_divr) { if (divr[i] < threshold_divr) {
status[i] |= FixRHEO::STATUS_SURFACE; status[i] |= STATUS_SURFACE;
rsurface[i] = 0.0; rsurface[i] = 0.0;
if (coordination[i] < threshold_z) if (coordination[i] < threshold_z)
status[i] |= FixRHEO::STATUS_SPLASH; status[i] |= STATUS_SPLASH;
} }
} }
} }
} else { } else {
for (i = 0; i < nall; i++) { for (i = 0; i < nall; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
status[i] |= FixRHEO::STATUS_BULK; status[i] |= STATUS_BULK;
rsurface[i] = cut; rsurface[i] = cut;
if (coordination[i] < divR_limit) { if (coordination[i] < threshold_divr) {
status[i] |= FixRHEO::STATUS_SURFACE; status[i] |= STATUS_SURFACE;
rsurface[i] = 0.0; rsurface[i] = 0.0;
if (coordination[i] < threshold_z) if (coordination[i] < threshold_z)
status[i] |= FixRHEO::STATUS_SPLASH; status[i] |= STATUS_SPLASH;
} }
} }
} }
@ -297,23 +298,23 @@ void ComputeRHEOSurface::compute_peratom()
dx[0] = xtmp - x[j][0]; dx[0] = xtmp - x[j][0];
dx[1] = ytmp - x[j][1]; dx[1] = ytmp - x[j][1];
dx[2] = ztmp - x[j][2]; dx[2] = ztmp - x[j][2];
rsq = lensq(dx); rsq = lensq3(dx);
if (rsq < cutsq) { if (rsq < cutsq) {
if ((status[i] & FixRHEO::STATUS_BULK) && (status[j] & FixRHEO::STATUS_SURFACE)) { if ((status[i] & STATUS_BULK) && (status[j] & STATUS_SURFACE)) {
status[i] &= FixRHEO::surfacemask; status[i] &= SURFACEMASK;
status[i] |= FixRHEO::STATUS_LAYER; status[i] |= STATUS_LAYER;
} }
if (status[j] & FixRHEO::STATUS_SURFACE) rsurface[i] = MIN(rsurface[i], sqrt(rsq)); if (status[j] & STATUS_SURFACE) rsurface[i] = MIN(rsurface[i], sqrt(rsq));
if (j < nlocal || newton) { if (j < nlocal || newton) {
if ((status[j] & FixRHEO::STATUS_BULK) && (status[i] & FixRHEO::STATUS_SURFACE)) { if ((status[j] & STATUS_BULK) && (status[i] & STATUS_SURFACE)) {
status[j] &= FixRHEO::surfacemask; status[j] &= SURFACEMASK;
status[j] |= FixRHEO::STATUS_LAYER; status[j] |= STATUS_LAYER;
} }
if (status[i] & FixRHEO::STATUS_SURFACE) rsurface[j] = MIN(rsurface[j], sqrt(rsq)); if (status[i] & STATUS_SURFACE) rsurface[j] = MIN(rsurface[j], sqrt(rsq));
} }
} }
} }
@ -371,7 +372,7 @@ void ComputeRHEOSurface::unpack_reverse_comm(int n, int *list, double *buf)
} else if (comm_stage == 1) { } else if (comm_stage == 1) {
temp = (int) buf[m++]; temp = (int) buf[m++];
if ((status[j] & FixRHEO::STATUS_BULK) && (temp & FixRHEO::STATUS_LAYER)) if ((status[j] & STATUS_BULK) && (temp & STATUS_LAYER))
status[j] = temp; status[j] = temp;
rsurface[j] = MIN(rsurface[j], buf[m++]); rsurface[j] = MIN(rsurface[j], buf[m++]);

10
src/RHEO/compute_rheo_surface.h
View File

@ -17,8 +17,8 @@ ComputeStyle(RHEO/SURFACE,ComputeRHEOSurface)
// clang-format on // clang-format on
#else #else
#ifndef LMP_COMPUTE_RHEO_INTERFACE_H #ifndef LMP_COMPUTE_RHEO_SURFACE_H
#define LMP_COMPUTE_RHEO_INTERFACE_H #define LMP_COMPUTE_RHEO_SURFACE_H
#include "compute.h" #include "compute.h"
@ -36,18 +36,18 @@ class ComputeRHEOSurface : public Compute {
int pack_forward_comm(int, int *, double *, int, int *) override; int pack_forward_comm(int, int *, double *, int, int *) override;
void unpack_forward_comm(int, int, double *) override; void unpack_forward_comm(int, int, double *) override;
double **nsurface, **rsurface; double **nsurface, *rsurface;
class FixRHEO *fix_rheo; class FixRHEO *fix_rheo;
private: private:
double cut, cutsq, rho0, threshold_divr; double cut, cutsq, rho0, threshold_divr;
int surface_style, nmax_old, threshold_z; int surface_style, nmax_store, threshold_z;
double **B, **gradC, *divr; double **B, **gradC, *divr;
int threshold_style, comm_stage; int threshold_style, comm_stage;
class NeighList *list; class NeighList *list;
class ComputeRHEOKernel *compute_kernel; class ComputeRHEOKernel *compute_kernel;
class ComputeRHEOSolids *compute_solids; class ComputeRHEOInterface *compute_interface;
}; };
} // namespace LAMMPS_NS } // namespace LAMMPS_NS

38
src/RHEO/compute_rheo_vshift.cpp
View File

@ -32,6 +32,7 @@
#include "neigh_request.h" #include "neigh_request.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace RHEO_NS;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -47,15 +48,15 @@ ComputeRHEOVShift::ComputeRHEOVShift(LAMMPS *lmp, int narg, char **arg) :
// Create vshift array if it doesn't already exist // Create vshift array if it doesn't already exist
// Create a custom atom property so it works with compute property/atom // 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 // Do not create grow callback as there's no reason to copy/exchange data
// Manually grow if nmax_old exceeded // Manually grow if nmax_store exceeded
int tmp1, tmp2; int tmp1, tmp2;
int index = atom->find_custom("rheo_vshift", tmp1, tmp2); int index = atom->find_custom("rheo_vshift", tmp1, tmp2);
if (index == -1) { if (index == -1) {
index = atom->add_custom("rheo_vshift", 1, 3); index = atom->add_custom("rheo_vshift", 1, 3);
nmax_old = atom->nmax; nmax_store = atom->nmax;
} }
vshift = atom->dvector[index]; vshift = atom->darray[index];
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -108,16 +109,17 @@ void ComputeRHEOVShift::compute_peratom()
int *jlist; int *jlist;
int inum, *ilist, *numneigh, **firstneigh; int inum, *ilist, *numneigh, **firstneigh;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
double **x = atom->x;
double **v = atom->v;
int *type = atom->type; int *type = atom->type;
int *status = atom->status; int *status = atom->status;
int *surface = atom->surface; int *mask = atom->mask;
double **x = atom->x;
double **v = atom->v;
double *rho = atom->rho; double *rho = atom->rho;
double *mass = atom->mass; double *mass = atom->mass;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
int newton_pair = force->newton_pair; int newton_pair = force->newton_pair;
inum = list->inum; inum = list->inum;
@ -125,9 +127,9 @@ void ComputeRHEOVShift::compute_peratom()
numneigh = list->numneigh; numneigh = list->numneigh;
firstneigh = list->firstneigh; firstneigh = list->firstneigh;
if (nmax_old < atom->nmax) { if (nmax_store < atom->nmax) {
memory->grow(vshift, atom->nmax, 3, "atom:rheo_vshift"); memory->grow(vshift, atom->nmax, 3, "atom:rheo_vshift");
nmax_old = atom->nmax; nmax_store = atom->nmax;
} }
for (i = 0; i < nall; i++) for (i = 0; i < nall; i++)
@ -143,15 +145,15 @@ void ComputeRHEOVShift::compute_peratom()
jlist = firstneigh[i]; jlist = firstneigh[i];
jnum = numneigh[i]; jnum = numneigh[i];
imass = mass[itype]; imass = mass[itype];
fluidi = status[i] & FixRHEO::STATUS_FLUID; fluidi = status[i] & STATUS_FLUID;
for (jj = 0; jj < jnum; jj++) { for (jj = 0; jj < jnum; jj++) {
j = jlist[jj]; j = jlist[jj];
j &= NEIGHMASK; j &= NEIGHMASK;
fluidj = status[j] & FixRHEO::STATUS_FLUID; fluidj = status[j] & STATUS_FLUID;
if ((!fluidi) && (!fluidj)) continue; if ((!fluidi) && (!fluidj)) continue;
if (!(status[i] & FixRHEO::STATUS_SHIFT) && !(status[j] & FixRHEO::STATUS_SHIFT)) continue; if (!(status[i] & STATUS_SHIFT) && !(status[j] & STATUS_SHIFT)) continue;
dx[0] = xtmp - x[j][0]; dx[0] = xtmp - x[j][0];
dx[1] = ytmp - x[j][1]; dx[1] = ytmp - x[j][1];
@ -175,10 +177,10 @@ void ComputeRHEOVShift::compute_peratom()
// Add corrections for walls // Add corrections for walls
if (fluidi && (!fluidj)) { if (fluidi && (!fluidj)) {
compute_interface->correct_v(v[i], v[j], vi, i, j); compute_interface->correct_v(vi, vj, i, j);
rhoj = compute_interface->correct_rho(j,i); rhoj = compute_interface->correct_rho(j,i);
} else if ((!fluidi) && fluidj) { } else if ((!fluidi) && fluidj) {
compute_interface->correct_v(v[j], v[i], vj, j, i); compute_interface->correct_v(vj, vi, j, i);
rhoi = compute_interface->correct_rho(i,j); rhoi = compute_interface->correct_rho(i,j);
} else if ((!fluidi) && (!fluidj)) { } else if ((!fluidi) && (!fluidj)) {
rhoi = 1.0; rhoi = 1.0;
@ -215,7 +217,7 @@ void ComputeRHEOVShift::compute_peratom()
} }
} }
if (newton_pair) comm->reverse_comm_compute(this); if (newton_pair) comm->reverse_comm(this);
} }
@ -239,7 +241,7 @@ void ComputeRHEOVShift::correct_surfaces()
double nx,ny,nz,vx,vy,vz; double nx,ny,nz,vx,vy,vz;
for (i = 0; i < nlocal; i++) { for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
if ((status[i] & FixRHEO::STATUS_SURFACE) || (status[i] & FixRHEO::STATUS_LAYER)) { if ((status[i] & STATUS_SURFACE) || (status[i] & STATUS_LAYER)) {
nx = nsurf[i][0]; nx = nsurf[i][0];
ny = nsurf[i][1]; ny = nsurf[i][1];
vx = vshift[i][0]; vx = vshift[i][0];
@ -297,6 +299,6 @@ void ComputeRHEOVShift::unpack_reverse_comm(int n, int *list, double *buf)
double ComputeRHEOVShift::memory_usage() double ComputeRHEOVShift::memory_usage()
{ {
double bytes = 3 * nmax_old * sizeof(double); double bytes = 3 * nmax_store * sizeof(double);
return bytes; return bytes;
} }

2
src/RHEO/compute_rheo_vshift.h
View File

@ -40,7 +40,7 @@ class ComputeRHEOVShift : public Compute {
class FixRHEO *fix_rheo; class FixRHEO *fix_rheo;
private: private:
int nmax_old; int nmax_store;
double dtv, cut, cutsq, cutthird; double dtv, cut, cutsq, cutthird;
int surface_flag; int surface_flag;

79
src/RHEO/fix_rheo.cpp
View File

@ -33,6 +33,7 @@
#include "utils.h" #include "utils.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace RHEO_NS;
using namespace FixConst; using namespace FixConst;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -68,6 +69,7 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR,"Insufficient arguments for fix rheo command"); error->all(FLERR,"Insufficient arguments for fix rheo command");
h = utils::numeric(FLERR,arg[3],false,lmp); h = utils::numeric(FLERR,arg[3],false,lmp);
cut = h;
if (strcmp(arg[4],"Quintic") == 0) { if (strcmp(arg[4],"Quintic") == 0) {
kernel_style = QUINTIC; kernel_style = QUINTIC;
} else if (strcmp(arg[4],"CRK0") == 0) { } else if (strcmp(arg[4],"CRK0") == 0) {
@ -101,7 +103,7 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg],"interface/reconstruction") == 0) { } else if (strcmp(arg[iarg],"interface/reconstruction") == 0) {
interface_flag = 1; interface_flag = 1;
} else if (strcmp(arg[iarg],"rhosum") == 0) { } else if (strcmp(arg[iarg],"rho/sum") == 0) {
rhosum_flag = 1; rhosum_flag = 1;
} else if (strcmp(arg[iarg],"rho0") == 0) { } else if (strcmp(arg[iarg],"rho0") == 0) {
if(iarg + 1 >= narg) error->all(FLERR,"Illegal rho0 option in fix rheo"); if(iarg + 1 >= narg) error->all(FLERR,"Illegal rho0 option in fix rheo");
@ -137,7 +139,8 @@ FixRHEO::~FixRHEO()
void FixRHEO::post_constructor() void FixRHEO::post_constructor()
{ {
compute_kernel = dynamic_cast<ComputeRHEOKernel *>(modify->add_compute("rheo_kernel all RHEO/KERNEL {}", kernel_style)); compute_kernel = dynamic_cast<ComputeRHEOKernel *>(modify->add_compute(
fmt::format("rheo_kernel all RHEO/KERNEL {}", kernel_style)));
compute_kernel->fix_rheo = this; 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";
@ -146,22 +149,26 @@ void FixRHEO::post_constructor()
compute_grad->fix_rheo = this; compute_grad->fix_rheo = this;
if (rhosum_flag) { if (rhosum_flag) {
compute_rhosum = dynamic_cast<ComputeRHEORhoSum *>(modify->add_compute("rheo_rho_sum all RHEO/RHO/SUM")); compute_rhosum = dynamic_cast<ComputeRHEORhoSum *>(modify->add_compute(
"rheo_rho_sum all RHEO/RHO/SUM"));
compute_rhosum->fix_rheo = this; compute_rhosum->fix_rheo = this;
} }
if (shift_flag) { 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; compute_vshift->fix_rheo = this;
} }
if (interface_flag) { 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(
"rheo_interface all RHEO/INTERFACE"));
compute_interface->fix_rheo = this; compute_interface->fix_rheo = this;
} }
if (surface_flag) { if (surface_flag) {
compute_surface = dynamic_cast<ComputeRHEOSurface *>(modify->add_compute(fmt::format("rheo_surface all RHEO/SURFACE"))); compute_surface = dynamic_cast<ComputeRHEOSurface *>(modify->add_compute(
"rheo_surface all RHEO/SURFACE"));
compute_surface->fix_rheo = this; compute_surface->fix_rheo = this;
} }
} }
@ -193,7 +200,7 @@ void FixRHEO::init()
void FixRHEO::setup_pre_force(int /*vflag*/) void FixRHEO::setup_pre_force(int /*vflag*/)
{ {
// Check to confirm accessory fixes do not preceed FixRHEO // Check to confirm accessory fixes do not preceed FixRHEO
// Note: these fixes set this flag in setup_pre_force() // Note: fixes set this flag in setup_pre_force()
if (viscosity_fix_defined || pressure_fix_defined || thermal_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"); error->all(FLERR, "Fix RHEO must be defined before all other RHEO fixes");
@ -206,23 +213,23 @@ void FixRHEO::setup_pre_force(int /*vflag*/)
void FixRHEO::setup(int /*vflag*/) void FixRHEO::setup(int /*vflag*/)
{ {
// Confirm all accessory fixes are defined // Confirm all accessory fixes are defined
// Note: these fixes set this flag in setup_pre_force() // Note: fixes set this flag in setup_pre_force()
if (!viscosity_fix_defined) if (!viscosity_fix_defined)
error->all(FLERR, "Missing fix rheo/viscosity"); error->all(FLERR, "Missing fix rheo/viscosity");
if (!pressure_fix_defined) if (!pressure_fix_defined)
error->all(FLERR, "Missing fix rheo/pressure"); error->all(FLERR, "Missing fix rheo/pressure");
if(!thermal_fix_defined && thermal_flag) if((!thermal_fix_defined) && thermal_flag)
error->all(FLERR, "Missing fix rheo/thermal"); error->all(FLERR, "Missing fix rheo/thermal");
// Reset to zero for next run // Reset to zero for future runs
thermal_fix_defined = 0; thermal_fix_defined = 0;
viscosity_fix_defined = 0; viscosity_fix_defined = 0;
pressure_fix_defined = 0; pressure_fix_defined = 0;
// Check fixes cover all atoms (doesnt ensure user covers atoms created midrun) // Check fixes cover all atoms (may still fail if atoms are created)
// (pressure is currently required to be group all) // FixRHEOPressure currently requires group all
auto visc_fixes = modify->get_fix_by_style("rheo/viscosity"); auto visc_fixes = modify->get_fix_by_style("rheo/viscosity");
auto therm_fixes = modify->get_fix_by_style("rheo/thermal"); auto therm_fixes = modify->get_fix_by_style("rheo/thermal");
@ -232,12 +239,12 @@ void FixRHEO::setup(int /*vflag*/)
int covered; int covered;
for (int i = 0; i < atom->nlocal; i++) { for (int i = 0; i < atom->nlocal; i++) {
covered = 0; covered = 0;
for (auto fix in visc_fixes) for (auto fix : visc_fixes)
if (mask[i] & fix->groupbit) covered = 1; if (mask[i] & fix->groupbit) covered = 1;
if (!covered) v_coverage_flag = 0; if (!covered) v_coverage_flag = 0;
if (thermal_flag) { if (thermal_flag) {
covered = 0; covered = 0;
for (auto fix in therm_fixes) for (auto fix : therm_fixes)
if (mask[i] & fix->groupbit) covered = 1; if (mask[i] & fix->groupbit) covered = 1;
if (!covered) v_coverage_flag = 0; if (!covered) v_coverage_flag = 0;
} }
@ -253,11 +260,12 @@ void FixRHEO::setup(int /*vflag*/)
void FixRHEO::initial_integrate(int /*vflag*/) void FixRHEO::initial_integrate(int /*vflag*/)
{ {
// update v and x and rho of atoms in group // update v, x and rho of atoms in group
int i, a, b; int i, a, b;
double dtfm, divu; double dtfm, divu;
int dim = domain->dimension;
int *type = atom->type;
int *mask = atom->mask;
int *status = atom->status; int *status = atom->status;
double **x = atom->x; double **x = atom->x;
double **v = atom->v; double **v = atom->v;
@ -266,16 +274,14 @@ void FixRHEO::initial_integrate(int /*vflag*/)
double *drho = atom->drho; double *drho = atom->drho;
double *mass = atom->mass; double *mass = atom->mass;
double *rmass = atom->rmass; double *rmass = atom->rmass;
int rmass_flag = atom->rmass_flag;
double **gradr = compute_grad->gradr; double **gradr = compute_grad->gradr;
double **gradv = compute_grad->gradv; double **gradv = compute_grad->gradv;
double **vshift; double **vshift;
if (shift_flag) compute_vshift->vshift; if (shift_flag) compute_vshift->vshift;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int rmass_flag = atom->rmass_flag;
int dim = domain->dimension;
if (igroup == atom->firstgroup) if (igroup == atom->firstgroup)
nlocal = atom->nfirst; nlocal = atom->nfirst;
@ -333,7 +339,7 @@ void FixRHEO::initial_integrate(int /*vflag*/)
// Shifting atoms // Shifting atoms
if (shift_flag) { if (shift_flag) {
compute_vshift->correct_surfaces(); // COuld this be moved to preforce after the surface fix runs? compute_vshift->correct_surfaces(); // Could this be moved to preforce after the surface fix runs?
for (i = 0; i < nlocal; i++) { for (i = 0; i < nlocal; i++) {
if (!(status[i] & STATUS_SHIFT)) continue; if (!(status[i] & STATUS_SHIFT)) continue;
@ -376,6 +382,7 @@ void FixRHEO::pre_force(int /*vflag*/)
compute_vshift->compute_peratom(); compute_vshift->compute_peratom();
// Remove extra shifting/no force options // Remove extra shifting/no force options
int *mask = atom->mask;
int *status = atom->status; int *status = atom->status;
int nall = atom->nlocal + atom->nghost; int nall = atom->nlocal + atom->nghost;
for (int i = 0; i < nall; i++) { for (int i = 0; i < nall; i++) {
@ -393,26 +400,28 @@ void FixRHEO::pre_force(int /*vflag*/)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixRHEO::final_integrate() { void FixRHEO::final_integrate()
int *status = atom->status; {
double **gradv = compute_grad->gradv;
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *rho = atom->rho;
double *drho = atom->drho;
int *type = atom->type;
int *mask = atom->mask;
double *mass = atom->mass;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) if (igroup == atom->firstgroup)
nlocal = atom->nfirst; nlocal = atom->nfirst;
double dtfm, divu; double dtfm, divu;
double *rmass = atom->rmass;
int rmass_flag = atom->rmass_flag;
int i, a; int i, a;
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double **gradv = compute_grad->gradv;
double *rho = atom->rho;
double *drho = atom->drho;
double *mass = atom->mass;
double *rmass = atom->rmass;
int *type = atom->type;
int *mask = atom->mask;
int *status = atom->status;
int rmass_flag = atom->rmass_flag;
int dim = domain->dimension; int dim = domain->dimension;
// Update velocity // Update velocity

51
src/RHEO/fix_rheo.h
View File

@ -39,33 +39,10 @@ class FixRHEO : public Fix {
void reset_dt() override; void reset_dt() override;
// Model parameters // Model parameters
double h, rho0, csq; double h, cut, rho0, csq;
int zmin_kernel, zmin_surface; int zmin_kernel, zmin_surface;
int kernel_style, surface_style; int kernel_style, surface_style;
double divr_surface; double divr_surface;
enum {QUINTIC, CRK0, CRK1, CRK2};
enum {COORDINATION, DIVR};
// Status variables
enum {
// Phase status
STATUS_FLUID = 1 << 0,
STATUS_REACTIVE = 1 << 1,
STATUS_SOLID = 1 << 2,
STATUS_FREEZING = 1 << 3,
// Surface status
STATUS_BULK = 1 << 4,
STATUS_LAYER = 1 << 5,
STATUS_SURFACE = 1 << 6,
STATUS_SPLASH = 1 << 7,
// Temporary status options - reset in preforce
STATUS_SHIFT = 1 << 8,
STATUS_NO_FORCE = 1 << 9
};
int phasemask = 0xFFFFFFF0;
int surfacemask = 0xFFFFFF0F;
// Accessory fixes/computes // Accessory fixes/computes
int thermal_flag; int thermal_flag;
@ -89,6 +66,32 @@ class FixRHEO : public Fix {
double dtv, dtf; double dtv, dtf;
}; };
namespace RHEO_NS {
enum {QUINTIC, CRK0, CRK1, CRK2};
enum {COORDINATION, DIVR};
// Status variables
enum Status{
// Phase status
STATUS_FLUID = 1 << 0,
STATUS_REACTIVE = 1 << 1,
STATUS_SOLID = 1 << 2,
STATUS_FREEZING = 1 << 3,
// Surface status
STATUS_BULK = 1 << 4,
STATUS_LAYER = 1 << 5,
STATUS_SURFACE = 1 << 6,
STATUS_SPLASH = 1 << 7,
// Temporary status options - reset in preforce
STATUS_SHIFT = 1 << 8,
STATUS_NO_FORCE = 1 << 9
};
int PHASEMASK = 0xFFFFFFF0;
int SURFACEMASK = 0xFFFFFF0F;
} // namespace RHEO_NS
} // namespace LAMMPS_NS } // namespace LAMMPS_NS
#endif #endif

25
src/RHEO/fix_rheo_pressure.cpp
View File

@ -33,6 +33,8 @@ using namespace LAMMPS_NS;
using namespace FixConst; using namespace FixConst;
enum {NONE, LINEAR, CUBIC, TAITWATER}; enum {NONE, LINEAR, CUBIC, TAITWATER};
static constexpr double SEVENTH = 1.0 / 7.0;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) : FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) :
@ -43,7 +45,7 @@ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) :
pressure_style = NONE; pressure_style = NONE;
comm_forward = 1; comm_forward = 1;
nmax_old = 0; nmax_store = 0;
// Currently can only have one instance of fix rheo/pressure // Currently can only have one instance of fix rheo/pressure
if (igroup != 0) if (igroup != 0)
@ -112,13 +114,13 @@ void FixRHEOPressure::setup_pre_force(int /*vflag*/)
// Create pressure array if it doesn't already exist // Create pressure array if it doesn't already exist
// Create a custom atom property so it works with compute property/atom // 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 // Do not create grow callback as there's no reason to copy/exchange data
// Manually grow if nmax_old exceeded // Manually grow if nmax_store exceeded
int tmp1, tmp2; int tmp1, tmp2;
int index = atom->find_custom("rheo_pressure", tmp1, tmp2); int index = atom->find_custom("rheo_pressure", tmp1, tmp2);
if (index == -1) { if (index == -1) {
index = atom->add_custom("rheo_pressure", 1, 0); index = atom->add_custom("rheo_pressure", 1, 0);
nmax_old = atom->nmax; nmax_store = atom->nmax;
} }
pressure = atom->dvector[index]; pressure = atom->dvector[index];
@ -139,13 +141,11 @@ void FixRHEOPressure::pre_force(int /*vflag*/)
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
if (nmax_old < atom->nmax) { if (nmax_store < atom->nmax) {
memory->grow(pressure, atom->nmax, "atom:rheo_pressure"); memory->grow(pressure, atom->nmax, "atom:rheo_pressure");
nmax_old = atom->nmax; nmax_store = atom->nmax;
} }
if (pressure_style == TAITWATER) inv7 = 1.0 / 7.0;
for (i = 0; i < nlocal; i++) { for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
if (pressure_style == LINEAR) { if (pressure_style == LINEAR) {
@ -156,7 +156,7 @@ void FixRHEOPressure::pre_force(int /*vflag*/)
} else if (pressure_style == TAITWATER) { } else if (pressure_style == TAITWATER) {
rho_ratio = rho[i] / rho0inv; rho_ratio = rho[i] / rho0inv;
rr3 = rho_ratio * rho_ratio * rho_ratio; rr3 = rho_ratio * rho_ratio * rho_ratio;
pressure[i] = csq * rho0 * inv7 * (rr3 * rr3 * rho_ratio - 1.0); pressure[i] = csq * rho0 * SEVENTH * (rr3 * rr3 * rho_ratio - 1.0);
} }
} }
} }
@ -194,9 +194,10 @@ void FixRHEOPressure::unpack_forward_comm(int n, int first, double *buf)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
double FixRHEOPressure::calculate_p(double rho) double FixRHEOPressure::calc_pressure(double rho)
{ {
double rho; double p, dr, rr3, rho_ratio;
if (pressure_style == LINEAR) { if (pressure_style == LINEAR) {
p = csq * (rho - rho0); p = csq * (rho - rho0);
} else if (pressure_style == CUBIC) { } else if (pressure_style == CUBIC) {
@ -205,7 +206,7 @@ double FixRHEOPressure::calculate_p(double rho)
} else if (pressure_style == TAITWATER) { } else if (pressure_style == TAITWATER) {
rho_ratio = rho / rho0inv; rho_ratio = rho / rho0inv;
rr3 = rho_ratio * rho_ratio * rho_ratio; rr3 = rho_ratio * rho_ratio * rho_ratio;
p = csq * rho0 * inv7 * (rr3 * rr3 * rho_ratio - 1.0); p = csq * rho0 * SEVENTH * (rr3 * rr3 * rho_ratio - 1.0);
} }
return rho; return rho;
} }
@ -215,6 +216,6 @@ double FixRHEOPressure::calculate_p(double rho)
double FixRHEOPressure::memory_usage() double FixRHEOPressure::memory_usage()
{ {
double bytes = 0.0; double bytes = 0.0;
bytes += (size_t) nmax_old * sizeof(double); bytes += (size_t) nmax_store * sizeof(double);
return bytes; return bytes;
} }

5
src/RHEO/fix_rheo_pressure.h
View File

@ -35,13 +35,14 @@ class FixRHEOPressure : public Fix {
int pack_forward_comm(int, int *, double *, int, int *) override; int pack_forward_comm(int, int *, double *, int, int *) override;
void unpack_forward_comm(int, int, double *) override; void unpack_forward_comm(int, int, double *) override;
double memory_usage() override; double memory_usage() override;
double calculate_p(double); double calc_pressure(double);
private: private:
double c_cubic, csq, rho0, rho0inv; double c_cubic, csq, rho0, rho0inv;
double *pressure; double *pressure;
int pressure_style; int pressure_style;
int first_flag, last_flag; int first_flag, last_flag;
int nmax_old; int nmax_store;
class FixRHEO *fix_rheo; class FixRHEO *fix_rheo;
}; };

57
src/RHEO/fix_rheo_thermal.cpp
View File

@ -22,6 +22,7 @@
#include "comm.h" #include "comm.h"
#include "compute_rheo_grad.h" #include "compute_rheo_grad.h"
#include "compute_rheo_vshift.h" #include "compute_rheo_vshift.h"
#include "domain.h"
#include "error.h" #include "error.h"
#include "fix_rheo.h" #include "fix_rheo.h"
#include "force.h" #include "force.h"
@ -31,14 +32,15 @@
#include "update.h" #include "update.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace RHEO_NS;
using namespace FixConst; using namespace FixConst;
enum {NONE, CONSTANT, TYPE}; enum {NONE, CONSTANT, TYPE};
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) : FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), Tc_type(nullptr), kappa_type(nullptr), cv_type(nullptr), Fix(lmp, narg, arg), fix_rheo(nullptr), compute_grad(nullptr), compute_vshift(nullptr),
conductivity(nullptr) Tc_type(nullptr), kappa_type(nullptr), cv_type(nullptr), conductivity(nullptr)
{ {
if (narg < 4) error->all(FLERR,"Illegal fix command"); if (narg < 4) error->all(FLERR,"Illegal fix command");
@ -47,7 +49,7 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
conductivity_style = NONE; conductivity_style = NONE;
comm_forward = 1; comm_forward = 1;
nmax_old = 0; nmax_store = 0;
int ntypes = atom->ntypes; int ntypes = atom->ntypes;
int iarg = 3; int iarg = 3;
@ -181,13 +183,13 @@ void FixRHEOThermal::setup_pre_force(int /*vflag*/)
// Identify whether this is the first/last instance of fix thermal // Identify whether this is the first/last instance of fix thermal
// First will grow arrays, last will communicate // First will grow arrays, last will communicate
first_flag = 0 first_flag = 0;
last_flag = 0; last_flag = 0;
int i = 0; int i = 0;
auto fixlist = modify->get_fix_by_style("rheo/thermal"); auto fixlist = modify->get_fix_by_style("rheo/thermal");
for (const auto &ifix : fixlist) { for (const auto &fix : fixlist) {
if (strcmp(ifix->id, id) == 0) break; if (strcmp(fix->id, id) == 0) break;
i++; i++;
} }
@ -197,13 +199,13 @@ void FixRHEOThermal::setup_pre_force(int /*vflag*/)
// Create conductivity array if it doesn't already exist // Create conductivity array if it doesn't already exist
// Create a custom atom property so it works with compute property/atom // 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 // Do not create grow callback as there's no reason to copy/exchange data
// Manually grow if nmax_old exceeded // Manually grow if nmax_store exceeded
int tmp1, tmp2; int tmp1, tmp2;
index = atom->find_custom("rheo_conductivity", tmp1, tmp2); int index = atom->find_custom("rheo_conductivity", tmp1, tmp2);
if (index== -1) { if (index== -1) {
index = atom->add_custom("rheo_conductivity", 1, 0); index = atom->add_custom("rheo_conductivity", 1, 0);
nmax_old = atom->nmax; nmax_store = atom->nmax;
} }
conductivity = atom->dvector[index]; conductivity = atom->dvector[index];
@ -217,13 +219,16 @@ void FixRHEOThermal::initial_integrate(int /*vflag*/)
{ {
// update temperature from shifting // update temperature from shifting
if (!fix_rheo->shift_flag) return; if (!fix_rheo->shift_flag) return;
int i; int i, a;
int *status = atom->status; int *status = atom->status;
int *mask = atom->mask;
double *temperature = atom->temperature;
double **gradt = compute_grad->gradt; double **gradt = compute_grad->gradt;
double **vshift = compute_vshift->array_atom; double **vshift = compute_vshift->array_atom;
int *mask = atom->mask;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int dim = domain->dimension;
if (igroup == atom->firstgroup) if (igroup == atom->firstgroup)
nlocal = atom->nfirst; nlocal = atom->nfirst;
@ -248,14 +253,14 @@ void FixRHEOThermal::post_integrate()
double *heatflow = atom->heatflow; double *heatflow = atom->heatflow;
double *rho = atom->rho; double *rho = atom->rho;
int *mask = atom->mask; int *mask = atom->mask;
int *type = aotm->type; int *type = atom->type;
double cvi, Tci, Ti; double cvi, Tci, Ti;
//Integrate temperature and check status //Integrate temperature and check status
for (int i = 0; i < atom->nlocal; i++) { for (int i = 0; i < atom->nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
if (status[i] == FixRHEO::FLUID_NO_FORCE) continue; if (status[i] & STATUS_NO_FORCE) continue;
cvi = calc_cv(i); cvi = calc_cv(i);
temperature[i] += dtf * heatflow[i] / cvi; temperature[i] += dtf * heatflow[i] / cvi;
@ -265,15 +270,15 @@ void FixRHEOThermal::post_integrate()
if (Tc_style == CONSTANT) { if (Tc_style == CONSTANT) {
Tci = Tc; Tci = Tc;
} else if (Tc_style == TYPE) { } else if (Tc_style == TYPE) {
Tci = Tc_type[type[i]]); Tci = Tc_type[type[i]];
} }
if (Ti > Tci) { if (Ti > Tci) {
status[i] &= FixRHEO::phasemask; status[i] &= PHASEMASK;
status[i] |= FixRHEO::STATUS_FLUID; status[i] |= STATUS_FLUID;
} else if (!(status[i] & FixRHEO::STATUS_SOLID)) } else if (!(status[i] & STATUS_SOLID)) {
status[i] &= FixRHEO::phasemask; status[i] &= PHASEMASK;
status[i] |= FixRHEO::STATUS_FREEZING; status[i] |= STATUS_FREEZING;
} }
} }
} }
@ -288,14 +293,13 @@ void FixRHEOThermal::post_neighbor()
{ {
int i; int i;
int *type = atom->type; int *type = atom->type;
double *conductivity = atom->dvector[index_cond];
int *mask = atom->mask; int *mask = atom->mask;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost; int nall = nlocal + atom->nghost;
if (first_flag && (nmax_old < atom->nmax)) { if (first_flag && (nmax_store < atom->nmax)) {
memory->grow(conductivity, atom->nmax, "atom:rheo_conductivity"); memory->grow(conductivity, atom->nmax, "atom:rheo_conductivity");
nmax_old = atom->nmax; nmax_store = atom->nmax;
} }
if (conductivity_style == CONSTANT) { if (conductivity_style == CONSTANT) {
@ -304,7 +308,6 @@ void FixRHEOThermal::post_neighbor()
} else if (conductivity_style == TYPE) { } else if (conductivity_style == TYPE) {
for (i = 0; i < nall; i++) for (i = 0; i < nall; i++)
if (mask[i] & groupbit) conductivity[i] = kappa_type[type[i]]; if (mask[i] & groupbit) conductivity[i] = kappa_type[type[i]];
}
} }
} }
@ -329,9 +332,9 @@ void FixRHEOThermal::pre_force(int /*vflag*/)
//int *mask = atom->mask; //int *mask = atom->mask;
//int nlocal = atom->nlocal; //int nlocal = atom->nlocal;
//if (first_flag && (nmax_old < atom->nmax)) { //if (first_flag && (nmax_store < atom->nmax)) {
// memory->grow(conductivity, atom->nmax, "atom:rheo_conductivity"); // memory->grow(conductivity, atom->nmax, "atom:rheo_conductivity");
// nmax_old = atom->nmax; // nmax_store = atom->nmax;
//} //}
//if (conductivity_style == TBD) { //if (conductivity_style == TBD) {
@ -358,7 +361,7 @@ void FixRHEOThermal::final_integrate()
//Integrate temperature and check status //Integrate temperature and check status
for (int i = 0; i < atom->nlocal; i++) { for (int i = 0; i < atom->nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
if (status[i] & FixRHEO::STATUS_NO_FORCE) continue; if (status[i] & STATUS_NO_FORCE) continue;
cvi = calc_cv(i); cvi = calc_cv(i);
temperature[i] += dtf * heatflow[i] / cvi; temperature[i] += dtf * heatflow[i] / cvi;
@ -447,6 +450,6 @@ void FixRHEOThermal::unpack_reverse_comm(int n, int *list, double *buf)
double FixRHEOThermal::memory_usage() double FixRHEOThermal::memory_usage()
{ {
double bytes = 0.0; double bytes = 0.0;
bytes += (size_t) nmax_old * sizeof(double); bytes += (size_t) nmax_store * sizeof(double);
return bytes; return bytes;
} }

6
src/RHEO/fix_rheo_thermal.h
View File

@ -31,7 +31,7 @@ class FixRHEOThermal : public Fix {
int setmask() override; int setmask() override;
void init() override; void init() override;
void setup_pre_force(int) override; void setup_pre_force(int) override;
void initial_integrate() override; void initial_integrate(int) override;
void post_integrate() override; void post_integrate() override;
void post_neighbor() override; void post_neighbor() override;
void pre_force(int) override; void pre_force(int) override;
@ -53,9 +53,11 @@ class FixRHEOThermal : public Fix {
int cv_style; int cv_style;
int conductivity_style; int conductivity_style;
int first_flag, last_flag; int first_flag, last_flag;
int nmax_old; int nmax_store;
class FixRHEO *fix_rheo; class FixRHEO *fix_rheo;
class ComputeRHEOGrad *compute_grad;
class ComputeRHEOVShift *compute_vshift;
double calc_cv(int); double calc_cv(int);
}; };

27
src/RHEO/fix_rheo_viscosity.cpp
View File

@ -44,7 +44,7 @@ FixRHEOViscosity::FixRHEOViscosity(LAMMPS *lmp, int narg, char **arg) :
viscosity_style = NONE; viscosity_style = NONE;
comm_forward = 0; comm_forward = 0;
nmax_old = 0; nmax_store = 0;
int ntypes = atom->ntypes; int ntypes = atom->ntypes;
int iarg = 3; int iarg = 3;
@ -89,7 +89,7 @@ FixRHEOViscosity::~FixRHEOViscosity()
// Remove custom property if it exists // Remove custom property if it exists
int tmp1, tmp2, index; int tmp1, tmp2, index;
index = atom->find_custom("rheo_viscosity", tmp1, tmp2); index = atom->find_custom("rheo_viscosity", tmp1, tmp2);
if (index != -1) atom->remove_custom(index_visc, 1, 0); if (index != -1) atom->remove_custom(index, 1, 0);
memory->destroy(eta_type); memory->destroy(eta_type);
} }
@ -123,13 +123,13 @@ void FixRHEOViscosity::setup_pre_force(int /*vflag*/)
// Identify whether this is the first/last instance of fix viscosity // Identify whether this is the first/last instance of fix viscosity
// First will grow arrays, last will communicate // First will grow arrays, last will communicate
first_flag = 0 first_flag = 0;
last_flag = 0; last_flag = 0;
int i = 0; int i = 0;
auto fixlist = modify->get_fix_by_style("rheo/viscosity"); auto fixlist = modify->get_fix_by_style("rheo/viscosity");
for (const auto &ifix : fixlist) { for (const auto &fix : fixlist) {
if (strcmp(ifix->id, id) == 0) break; if (strcmp(fix->id, id) == 0) break;
i++; i++;
} }
@ -139,13 +139,13 @@ void FixRHEOViscosity::setup_pre_force(int /*vflag*/)
// Create viscosity array if it doesn't already exist // Create viscosity array if it doesn't already exist
// Create a custom atom property so it works with compute property/atom // 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 // Do not create grow callback as there's no reason to copy/exchange data
// Manually grow if nmax_old exceeded // Manually grow if nmax_store exceeded
int tmp1, tmp2; int tmp1, tmp2;
int index = atom->find_custom("rheo_viscosity", tmp1, tmp2); int index = atom->find_custom("rheo_viscosity", tmp1, tmp2);
if (index_visc == -1) { if (index == -1) {
index = atom->add_custom("rheo_viscosity", 1, 0); index = atom->add_custom("rheo_viscosity", 1, 0);
nmax_old = atom->nmax; nmax_store = atom->nmax;
} }
viscosity = atom->dvector[index]; viscosity = atom->dvector[index];
@ -167,9 +167,9 @@ void FixRHEOViscosity::post_neighbor()
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost; int nall = nlocal + atom->nghost;
if (first_flag && (nmax_old < atom->nmax)) { if (first_flag && (nmax_store < atom->nmax)) {
memory->grow(viscosity, atom->nmax, "atom:rheo_viscosity"); memory->grow(viscosity, atom->nmax, "atom:rheo_viscosity");
nmax_old = atom->nmax; nmax_store = atom->nmax;
} }
if (viscosity_style == CONSTANT) { if (viscosity_style == CONSTANT) {
@ -178,7 +178,6 @@ void FixRHEOViscosity::post_neighbor()
} else if (viscosity_style == TYPE) { } else if (viscosity_style == TYPE) {
for (i = 0; i < nall; i++) for (i = 0; i < nall; i++)
if (mask[i] & groupbit) viscosity[i] = eta_type[type[i]]; if (mask[i] & groupbit) viscosity[i] = eta_type[type[i]];
}
} }
} }
@ -197,9 +196,9 @@ void FixRHEOViscosity::pre_force(int /*vflag*/)
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int dim = domain->dimension; int dim = domain->dimension;
if (first_flag && (nmax_old < atom->nmax)) { if (first_flag && (nmax_store < atom->nmax)) {
memory->grow(viscosity, atom->nmax, "atom:rheo_viscosity"); memory->grow(viscosity, atom->nmax, "atom:rheo_viscosity");
nmax_old = atom->nmax; nmax_store = atom->nmax;
} }
if (viscosity_style == POWER) { if (viscosity_style == POWER) {
@ -260,6 +259,6 @@ void FixRHEOViscosity::unpack_forward_comm(int n, int first, double *buf)
double FixRHEOViscosity::memory_usage() double FixRHEOViscosity::memory_usage()
{ {
double bytes = 0.0; double bytes = 0.0;
bytes += (size_t) nmax_old * sizeof(double); bytes += (size_t) nmax_store * sizeof(double);
return bytes; return bytes;
} }

2
src/RHEO/fix_rheo_viscosity.h
View File

@ -43,7 +43,7 @@ class FixRHEOViscosity : public Fix {
double *viscosity; double *viscosity;
int viscosity_style; int viscosity_style;
int first_flag, last_flag; int first_flag, last_flag;
int nmax_old; int nmax_store;
class FixRHEO *fix_rheo; class FixRHEO *fix_rheo;
class ComputeRHEOGrad *compute_grad; class ComputeRHEOGrad *compute_grad;

39
src/RHEO/pair_rheo.cpp
View File

@ -39,9 +39,10 @@
#include <cmath> #include <cmath>
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace RHEO_NS;
using namespace MathExtra; using namespace MathExtra;
#define EPSILON 1e-2 static constexpr double EPSILON = 1e-2;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -83,6 +84,10 @@ void PairRHEO::compute(int eflag, int vflag)
int *ilist, *jlist, *numneigh, **firstneigh; int *ilist, *jlist, *numneigh, **firstneigh;
double imass, jmass, rsq, r, rinv; double imass, jmass, rsq, r, rinv;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
int dim = domain->dimension;
ev_init(eflag, vflag); ev_init(eflag, vflag);
double **gradv = compute_grad->gradv; double **gradv = compute_grad->gradv;
@ -123,11 +128,6 @@ void PairRHEO::compute(int eflag, int vflag)
conductivity = atom->dvector[index]; conductivity = atom->dvector[index];
} }
int *ilist, *jlist, *numneigh, **firstneigh;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
int dim = domain->dimension;
inum = list->inum; inum = list->inum;
ilist = list->ilist; ilist = list->ilist;
numneigh = list->numneigh; numneigh = list->numneigh;
@ -145,7 +145,7 @@ void PairRHEO::compute(int eflag, int vflag)
jnum = numneigh[i]; jnum = numneigh[i];
imass = mass[itype]; imass = mass[itype];
etai = viscosity[i]; etai = viscosity[i];
fluidi = status[i] & FixRHEO::STATUS_FLUID; fluidi = status[i] & STATUS_FLUID;
if (thermal_flag) { if (thermal_flag) {
kappai = conductivity[i]; kappai = conductivity[i];
Ti = temperature[i]; Ti = temperature[i];
@ -167,7 +167,7 @@ void PairRHEO::compute(int eflag, int vflag)
jmass = mass[jtype]; jmass = mass[jtype];
etaj = viscosity[j]; etaj = viscosity[j];
fluidj = status[j] & FixRHEO::STATUS_FLUID; fluidj = status[j] & STATUS_FLUID;
if (thermal_flag) { if (thermal_flag) {
Tj = temperature[j]; Tj = temperature[j];
kappaj = conductivity[j]; kappaj = conductivity[j];
@ -202,7 +202,7 @@ void PairRHEO::compute(int eflag, int vflag)
if (fluidi && (!fluidj)) { if (fluidi && (!fluidj)) {
compute_interface->correct_v(vi, vj, i, j); compute_interface->correct_v(vi, vj, i, j);
rhoj = compute_interface->correct_rho(j, i); rhoj = compute_interface->correct_rho(j, i);
Pj = fix_pressure->calculate_p(rhoj); Pj = fix_pressure->calc_pressure(rhoj);
if ((chi[j] > 0.9) && (r < (h * 0.5))) if ((chi[j] > 0.9) && (r < (h * 0.5)))
fmag = (chi[j] - 0.9) * (h * 0.5 - r) * rho0 * csq * h * rinv; fmag = (chi[j] - 0.9) * (h * 0.5 - r) * rho0 * csq * h * rinv;
@ -210,9 +210,9 @@ void PairRHEO::compute(int eflag, int vflag)
} else if ((!fluidi) && fluidj) { } else if ((!fluidi) && fluidj) {
compute_interface->correct_v(vj, vi, j, i); compute_interface->correct_v(vj, vi, j, i);
rhoi = compute_interface->correct_rho(i, j); rhoi = compute_interface->correct_rho(i, j);
Pi = calc_pressure(rhoi, itype); Pi = fix_pressure->calc_pressure(rhoi);
if (chi[i] > 0.9 && r < (h * 0.5)) { if (chi[i] > 0.9 && r < (h * 0.5))
fmag = (chi[i] - 0.9) * (h * 0.5 - r) * rho0 * csq * h * rinv; fmag = (chi[i] - 0.9) * (h * 0.5 - r) * rho0 * csq * h * rinv;
} else if ((!fluidi) && (!fluidj)) { } else if ((!fluidi) && (!fluidj)) {
@ -244,7 +244,7 @@ void PairRHEO::compute(int eflag, int vflag)
//Hydrostatic pressure forces //Hydrostatic pressure forces
fp_prefactor = voli * volj * (Pj + Pi); fp_prefactor = voli * volj * (Pj + Pi);
sub3(v1, vj, du); sub3(vi, vj, du);
//Add artificial viscous pressure if required //Add artificial viscous pressure if required
if (artificial_visc_flag && pair_avisc_flag){ if (artificial_visc_flag && pair_avisc_flag){
@ -423,15 +423,11 @@ void PairRHEO::setup()
if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use pair rheo"); if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use pair rheo");
fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]); fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
// Currently only allow one instance of fix rheo/pressure
fixes = modify->get_fix_by_style("rheo/pressure"); fixes = modify->get_fix_by_style("rheo/pressure");
if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo/pressure to use pair rheo"); if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo/pressure to use pair rheo");
fix_pressure = dynamic_cast<FixRHEOPressure *>(fixes[0]); fix_pressure = dynamic_cast<FixRHEOPressure *>(fixes[0]);
int tmp1, tmp2;
index_pressure = atom->find_custom("rheo_pressure", tmp1, tmp2);
if (index_pressure == -1) index_pressure = atom->add_custom("rheo_pressure", 1, 0);
else error->all(FLERR, "Cannot find pressure value in pair rheo");
compute_kernel = fix_rheo->compute_kernel; compute_kernel = fix_rheo->compute_kernel;
compute_grad = fix_rheo->compute_grad; compute_grad = fix_rheo->compute_grad;
compute_interface = fix_rheo->compute_interface; compute_interface = fix_rheo->compute_interface;
@ -449,9 +445,9 @@ void PairRHEO::setup()
error->all(FLERR,"Pair RHEO requires ghost atoms store velocity"); error->all(FLERR,"Pair RHEO requires ghost atoms store velocity");
if (laplacian_order == -1) { if (laplacian_order == -1) {
if (fix_rheo->kernel_type == FixRHEO::CRK2) if (fix_rheo->kernel_style == CRK2)
laplacian_order = 2; laplacian_order = 2;
else if (fix_rheo->kernel_type == FixRHEO::CRK1) else if (fix_rheo->kernel_style == CRK1)
laplacian_order = 1; laplacian_order = 1;
else else
laplacian_order = 0; laplacian_order = 0;
@ -468,8 +464,5 @@ double PairRHEO::init_one(int i, int j)
error->all(FLERR,"All pair rheo coeffs are not set"); error->all(FLERR,"All pair rheo coeffs are not set");
} }
cut[i][j] = h; return h;
cut[j][i] = cut[i][j];
return cut[i][j];
} }