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

make compatible with per-type and per-atom masses

Browse Source
This commit is contained in:
Axel Kohlmeyer
2024-07-30 16:16:05 -04:00
parent 1254d579f7
commit 5016a0848f
7 changed files with 47 additions and 21 deletions
Download Patch File Download Diff File Expand all files Collapse all files

10
src/RHEO/compute_rheo_grad.cpp
View File

@ -147,6 +147,7 @@ void ComputeRHEOGrad::compute_peratom()
int *status = atom->rheo_status;
int *type = atom->type;
double *mass = atom->mass;
double *rmass = atom->rmass;
int newton = force->newton;
int dim = domain->dimension;
@ -225,8 +226,13 @@ void ComputeRHEOGrad::compute_peratom()
}
}
Voli = mass[itype] / rhoi;
Volj = mass[jtype] / rhoj;
if (rmass) {
Voli = rmass[i] / rhoi;
Volj = rmass[j] / rhoj;
} else {
Voli = mass[itype] / rhoi;
Volj = mass[jtype] / rhoj;
}
vij[0] = vi[0] - vj[0];
vij[1] = vi[1] - vj[1];

20
src/RHEO/compute_rheo_interface.cpp
View File

@ -335,6 +335,7 @@ void ComputeRHEOInterface::store_forces()
int *type = atom->type;
int *mask = atom->mask;
double *mass = atom->mass;
double *rmass = atom->rmass;
double **f = atom->f;
// When this is called, fp_store stores the pressure force
@ -346,7 +347,8 @@ void ComputeRHEOInterface::store_forces()
if (fixlist.size() != 0) {
for (const auto &fix : fixlist) {
for (int i = 0; i < atom->nlocal; i++) {
minv = 1.0 / mass[type[i]];
if (rmass) minv = 1.0 / rmass[i];
else minv = 1.0 / mass[type[i]];
if (mask[i] & fix->groupbit)
for (int a = 0; a < 3; a++)
fp_store[i][a] = f[i][a] * minv;
@ -356,10 +358,18 @@ void ComputeRHEOInterface::store_forces()
}
}
} else {
for (int i = 0; i < atom->nlocal; i++) {
minv = 1.0 / mass[type[i]];
for (int a = 0; a < 3; a++)
fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv;
if (rmass) {
for (int i = 0; i < atom->nlocal; i++) {
minv = 1.0 / rmass[i];
for (int a = 0; a < 3; a++)
fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv;
}
} else {
for (int i = 0; i < atom->nlocal; i++) {
minv = 1.0 / mass[type[i]];
for (int a = 0; a < 3; a++)
fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv;
}
}
}

10
src/RHEO/compute_rheo_kernel.cpp
View File

@ -48,8 +48,6 @@ using namespace RHEO_NS;
using namespace MathConst;
using namespace MathExtra;
static constexpr int DELTA = 2000;
/* ---------------------------------------------------------------------- */
ComputeRHEOKernel::ComputeRHEOKernel(LAMMPS *lmp, int narg, char **arg) :
@ -581,6 +579,7 @@ void ComputeRHEOKernel::compute_peratom()
double **x = atom->x;
int *type = atom->type;
double *mass = atom->mass;
double *rmass = atom->rmass;
double *rho = atom->rho;
int *status = atom->rheo_status;
tagint *tag = atom->tag;
@ -626,7 +625,8 @@ void ComputeRHEOKernel::compute_peratom()
if (status[j] & PHASECHECK)
rhoj = compute_interface->correct_rho(j, i);
vj = mass[type[j]] / rhoj;
if (rmass) vj = rmass[j] / rhoj;
else vj = mass[type[j]] / rhoj;
M += w * vj;
}
}
@ -647,7 +647,6 @@ void ComputeRHEOKernel::compute_peratom()
jlist = firstneigh[i];
jnum = numneigh[i];
itype = type[i];
// Zero upper-triangle M and cut (will be symmetric):
for (a = 0; a < Mdim; a++) {
@ -675,7 +674,8 @@ void ComputeRHEOKernel::compute_peratom()
if (status[j] & PHASECHECK)
rhoj = compute_interface->correct_rho(j, i);
vj = mass[type[j]] / rhoj;
if (rmass) vj = rmass[j] / rhoj;
else vj = mass[type[j]] / rhoj;
//Populate the H-vector of polynomials (2D)
if (dim == 2) {

11
src/RHEO/compute_rheo_surface.cpp
View File

@ -113,6 +113,7 @@ void ComputeRHEOSurface::compute_peratom()
int *mask = atom->mask;
int *type = atom->type;
double *mass = atom->mass;
double *rmass = atom->rmass;
double *rho = atom->rho;
int *coordination = compute_kernel->coordination;
@ -173,9 +174,13 @@ void ComputeRHEOSurface::compute_peratom()
}
}
Voli = mass[itype] / rhoi;
Volj = mass[jtype] / rhoj;
if (rmass) {
Voli = rmass[i] / rhoi;
Volj = rmass[j] / rhoj;
} else {
Voli = mass[itype] / rhoi;
Volj = mass[jtype] / rhoj;
}
compute_kernel->calc_dw_quintic(i, j, dx[0], dx[1], dx[2], sqrt(rsq), dWij, dWji);
for (a = 0; a < dim; a++) {

7
src/RHEO/compute_rheo_vshift.cpp
View File

@ -108,6 +108,7 @@ void ComputeRHEOVShift::compute_peratom()
double **v = atom->v;
double *rho = atom->rho;
double *mass = atom->mass;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
@ -140,7 +141,8 @@ void ComputeRHEOVShift::compute_peratom()
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
imass = mass[itype];
if (rmass) imass = rmass[i];
else imass = mass[itype];
fluidi = !(status[i] & PHASECHECK);
for (jj = 0; jj < jnum; jj++) {
@ -160,7 +162,8 @@ void ComputeRHEOVShift::compute_peratom()
if (rsq < cutsq) {
jtype = type[j];
jmass = mass[jtype];
if (rmass) jmass = rmass[j];
else jmass = mass[jtype];
r = sqrt(rsq);
rinv = 1 / r;

3
src/RHEO/fix_rheo.cpp
View File

@ -469,7 +469,6 @@ void FixRHEO::final_integrate()
int *mask = atom->mask;
int *status = atom->rheo_status;
int rmass_flag = atom->rmass_flag;
int dim = domain->dimension;
// Update velocity
@ -477,7 +476,7 @@ void FixRHEO::final_integrate()
if (mask[i] & groupbit) {
if (status[i] & STATUS_NO_INTEGRATION) continue;
if (rmass_flag) {
if (rmass) {
dtfm = dtf / rmass[i];
} else {
dtfm = dtf / mass[type[i]];

7
src/RHEO/pair_rheo.cpp
View File

@ -103,6 +103,7 @@ void PairRHEO::compute(int eflag, int vflag)
double **f = atom->f;
double *rho = atom->rho;
double *mass = atom->mass;
double *rmass = atom->rmass;
double *drho = atom->drho;
double *pressure = atom->pressure;
double *viscosity = atom->viscosity;
@ -149,7 +150,8 @@ void PairRHEO::compute(int eflag, int vflag)
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
imass = mass[itype];
if (rmass) imass = rmass[i];
else imass = mass[itype];
etai = viscosity[i];
fluidi = !(status[i] & PHASECHECK);
if (thermal_flag) {
@ -171,7 +173,8 @@ void PairRHEO::compute(int eflag, int vflag)
r = sqrt(rsq);
rinv = 1 / r;
jmass = mass[jtype];
if (rmass) jmass = rmass[i];
else jmass = mass[jtype];
etaj = viscosity[j];
fluidj = !(status[j] & PHASECHECK);
if (thermal_flag) {