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

Revising interface compute

Browse Source
This commit is contained in:
jtclemm
2023-04-15 21:15:37 -06:00
parent 986cfd6641
commit 5980fdf9fd
6 changed files with 126 additions and 120 deletions
Download Patch File Download Diff File Expand all files Collapse all files

185
src/RHEO/compute_rheo_interface.cpp
View File

@ -61,8 +61,8 @@ ComputeRHEOInterface::~ComputeRHEOInterface()
if (index != -1) atom->remove_custom(index, 1, 0);
if (id_fix_pa && modify->nfix) modify->delete_fix(id_fix_pa);
delete[] id_fix_pa;
memory->destroy(fx_m_norm);
memory->destroy(norm);
memory->destroy(normwf);
}
@ -72,7 +72,10 @@ ComputeRHEOInterface::~ComputeRHEOInterface()
void ComputeRHEOInterface::init()
{
compute_kernel = fix_rheo->compute_kernel;
rho0 = fix_rheo->rho0;
cut = fix_rheo->cut;
cs = fix_rheo->cs;
cs_inv = 1.0 / cs;
cutsq = cut * cut;
wall_max = sqrt(3.0) / 12.0 * cut;
@ -87,21 +90,25 @@ void ComputeRHEOInterface::init()
int index = atom->find_custom("rheo_chi", tmp1, tmp2);
if (index == -1) {
index = atom->add_custom("rheo_chi", 1, 0);
memory->destroy(norm);
memory->destroy(normwf);
memory->create(norm, nmax, "rheo/interface:norm");
memory->create(normwf, nmax, "rheo/interface:normwf");
nmax_old = nmax;
}
chi = atom->dvector[index];
// For fpressure, go ahead and create an instance of fix property atom
// For fp_store, go ahead and create an instance of fix property atom
// Need restarts + exchanging with neighbors since it needs to persist
// between timesteps (fix property atom will handle callbacks)
index = atom->find_custom("rheo_pressure", tmp1, tmp2);
index = atom->find_custom("fp_store", tmp1, tmp2);
if (index == -1) {
id_fix_pa = utils::strdup(id + std::string("_fix_property_atom"));
modify->add_fix(fmt::format("{} all property/atom d2_f_pressure 3", id_fix_pa)));
index = atom->find_custom("rheo_pressure", tmp1, tmp2);
modify->add_fix(fmt::format("{} all property/atom d2_fp_store 3", id_fix_pa)));
index = atom->find_custom("fp_store", tmp1, tmp2);
}
f_pressure = atom->darray[index];
fp_store = atom->darray[index];
// need an occasional half neighbor list
neighbor->add_request(this, NeighConst::REQ_HALF);
@ -116,51 +123,37 @@ void ComputeRHEOInterface::init_list(int /*id*/, NeighList *ptr)
/* ---------------------------------------------------------------------- */
// Left off here
void ComputeRHEOInterface::compute_peratom()
{
int i, j, ii, jj, jnum, itype, jtype, phase_match;
double xtmp, ytmp, ztmp, delx, dely, delz, rsq;
int *jlist;
double w;
int i, j, ii, jj, jnum, itype, jtype, fluidi, fluidj, status_match;
double xtmp, ytmp, ztmp, delx, dely, delz, rsq, w, dot;
// neighbor list ariables
int inum, *ilist, *numneigh, **firstneigh;
int inum, *ilist, *jlist, *numneigh, **firstneigh;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
double **x = atom->x;
int *type = atom->type;
int newton = force->newton;
int *phase = atom->phase;
int *status = atom->status;
double *rho = atom->rho;
double *fx = atom->dvector[index_fx];
double *fy = atom->dvector[index_fy];
double *fz = atom->dvector[index_fz];
double mi, mj;
//Declare mass pointer to calculate acceleration from force
double *mass = atom->mass;
cs2 = 1.0;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
if (atom->nmax > nmax) {
nmax = atom->nmax;
fix_chi->grow_arrays(nmax);
chi = fix_chi->vstore;
if (atom->nmax > nmax_old) {
nmax_old = atom->nmax;
memory->destroy(norm);
memory->destroy(normwf);
memory->create(norm, nmax, "RHEO/chi:norm");
memory->create(normwf, nmax, "RHEO/chi:normwf");
memory->create(norm, nmax_old, "rheo/interface:norm");
memory->create(normwf, nmax_old, "rheo/interface:normwf");
memory->grow(chi, nmax_old, "rheo/interface:chi");
}
for (i = 0; i < nall; i++) {
if (phase[i] > FixRHEO::FLUID_MAX) rho[i] = 0.0;
if (!(status[i] & FixRHEO::STATUS_FLUID)) rho[i] = 0.0;
normwf[i] = 0.0;
norm[i] = 0.0;
chi[i] = 0.0;
@ -172,9 +165,9 @@ void ComputeRHEOInterface::compute_peratom()
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
fluidi = status[i] & FixRHEO::STATUS_FLUID;
jlist = firstneigh[i];
jnum = numneigh[i];
mi = mass[type[i]];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
@ -187,37 +180,36 @@ void ComputeRHEOInterface::compute_peratom()
if (rsq < cutsq) {
jtype = type[j];
mj = mass[type[j]];
fluidj = status[j] & FixRHEO::STATUS_FLUID;
w = compute_kernel->calc_w_quintic(i, j, delx, dely, delz, sqrt(rsq));
phase_match = 0;
status_match = 0;
norm[i] += w;
if ((phase[i] <= FixRHEO::FLUID_MAX && phase[j] <= FixRHEO::FLUID_MAX)
|| (phase[i] > FixRHEO::FLUID_MAX && phase[j] > FixRHEO::FLUID_MAX)) {
phase_match = 1;
}
if ((fluidi && fluidj) || ((!fluid) && (!fluidj)))
status_match = 1;
if (phase_match) {
if (status_match) {
chi[i] += w;
} else {
if (phase[i] > FixRHEO::FLUID_MAX) {
//speed of sound and rho0 assumed to = 1 (units in density not pressure)
//In general, rho is calculated using the force vector on the wall particle
//fx stores f-fp
rho[i] += w*(cs2*(rho[j] - 1.0) - rho[j]*((-fx[j]/mj+fx[i]/mi)*delx + (-fy[j]/mj+fy[i]/mi)*dely + (-fz[j]/mj+fz[i]/mi)*delz));
//For the specific taste case whre force on wall particles = 0
//rho[i] += w*(1.0*(rho[j] - 1.0) + rho[j]*(1e-3*delx));
if (!fluidi) {
dot = (-fp_store[0][j] + fp_store[0][i]) * delx;
dot += (-fp_store[1][j] + fp_store[1][i]) * dely;
dot += (-fp_store[2][j] + fp_store[2][i]) * delz;
rho[i] += w * (cs * (rho[j] - rho0) - rho[j] * dot);
normwf[i] += w;
}
}
if (newton || j < nlocal) {
norm[j] += w;
if (phase_match) {
if (status_match) {
chi[j] += w;
} else {
if (phase[j] > FixRHEO::FLUID_MAX) {
rho[j] += w*(cs2*(rho[i] - 1.0) + rho[i]*((-fx[i]/mi+fx[j]/mj)*delx + (-fy[i]/mi+fy[j]/mj)*dely + (-fz[i]/mi+fz[j]/mj)*delz));
if (!fluidj) {
dot = (-fp_store[0][i] + fp_store[0][j]) * delx;
dot += (-fp_store[1][i] + fp_store[1][j]) * dely;
dot += (-fp_store[2][i] + fp_store[2][j]) * delz;
rho[j] += w * (cs * (rho[i] - rho0) + rho[i] * dot);
normwf[j] += w;
}
}
@ -226,41 +218,41 @@ void ComputeRHEOInterface::compute_peratom()
}
}
if (newton) comm->reverse_comm_compute(this);
if (newton) comm->reverse_comm(this);
for (i = 0; i < nlocal; i++) {
if (norm[i] != 0.0) chi[i] /= norm[i];
if (normwf[i] != 0.0) { // Only if it's a wall particle
rho[i] = 1.0 + (rho[i] / normwf[i])/cs2; // Stores rho for solid particles 1+Pw in Adami Adams 2012
if (rho[i] < EPSILON) rho[i] = EPSILON;
}
if (normwf[i] == 0.0 && phase[i] > FixRHEO::FLUID_MAX) rho[i] = 1.0;
// Recalculate rho for non-fluid particles
if (!(status[i] & FixRHEO::STATUS_FLUID)) {
if (normwf[i] != 0.0) {
// Stores rho for solid particles 1+Pw in Adami Adams 2012
rho[i] = MAX(EPSILON, rho0 + (rho[i] / normwf[i]) * cs_inv);
} else {
rho[i] = rho0;
}
}
}
comm_stage = 1;
comm_forward = 2;
comm->forward_comm_compute(this);
comm->forward_comm(this);
}
/* ---------------------------------------------------------------------- */
int ComputeRHEOInterface::pack_forward_comm(int n, int *list, double *buf,
int /*pbc_flag*/, int * /*pbc*/)
int ComputeRHEOInterface::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/)
{
int i,j,k,m;
m = 0;
double *rho = atom->rho;
double *fx = atom->dvector[index_fx];
double *fy = atom->dvector[index_fy];
double *fz = atom->dvector[index_fz];
for (i = 0; i < n; i++) {
j = list[i];
if (comm_stage == 0) {
buf[m++] = fx[j];
buf[m++] = fy[j];
buf[m++] = fz[j];
buf[m++] = fp_store[j][0];
buf[m++] = fp_store[j][1];
buf[m++] = fp_store[j][2];
} else {
buf[m++] = chi[j];
buf[m++] = rho[j];
@ -275,20 +267,16 @@ void ComputeRHEOInterface::unpack_forward_comm(int n, int first, double *buf)
{
int i, k, m, last;
double *rho = atom->rho;
double *fx = atom->dvector[index_fx];
double *fy = atom->dvector[index_fy];
double *fz = atom->dvector[index_fz];
m = 0;
last = first + n;
for (i = first; i < last; i++) {
if (comm_stage == 0) {
fx[i] = buf[m++];
fy[i] = buf[m++];
fz[i] = buf[m++];
fp_store[i][0] = buf[m++];
fp_store[i][1] = buf[m++];
fp_store[i][2] = buf[m++];
} else {
chi[i] = buf[m++];
rho[i] = buf[m++]; // Won't do anything for fluids
rho[i] = buf[m++];
}
}
}
@ -317,13 +305,13 @@ void ComputeRHEOInterface::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,k,j,m;
double *rho = atom->rho;
int *phase = atom->phase;
int *status = atom->status;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
norm[j] += buf[m++];
chi[j] += buf[m++];
if (phase[j] > FixRHEO::FLUID_MAX){
if (!(status[j] & FixRHEO::STATUS_FLUID)){
normwf[j] += buf[m++];
rho[j] += buf[m++];
} else {
@ -339,22 +327,23 @@ void ComputeRHEOInterface::correct_v(double *vi, double *vj, double *vi_out, int
{
double wall_prefactor, wall_denom, wall_numer;
wall_numer = 2.0*cut*(chi[i]-0.5);
wall_numer = 2.0 * cut * (chi[i] - 0.5);
if (wall_numer < 0) wall_numer = 0;
wall_denom = 2.0*cut*(chi[j]-0.5);
wall_denom = 2.0 * cut * (chi[j] - 0.5);
if (wall_denom < wall_max) wall_denom = wall_max;
wall_prefactor = wall_numer / wall_denom;
vi_out[0] = (vi[0]-vj[0])*wall_prefactor + vi[0];
vi_out[1] = (vi[1]-vj[1])*wall_prefactor + vi[1];
vi_out[2] = (vi[2]-vj[2])*wall_prefactor + vi[2];
vi_out[0] = (vi[0] - vj[0]) * wall_prefactor + vi[0];
vi_out[1] = (vi[1] - vj[1]) * wall_prefactor + vi[1];
vi_out[2] = (vi[2] - vj[2]) * wall_prefactor + vi[2];
}
/* ---------------------------------------------------------------------- */
double ComputeRHEOInterface::correct_rho(int i, int j) // i is wall, j is fluid
double ComputeRHEOInterface::correct_rho(int i, int j)
{
// i is wall, j is fluid
//In future may depend on atom type j's pressure equation
return atom->rho[i];
}
@ -363,42 +352,46 @@ double ComputeRHEOInterface::correct_rho(int i, int j) // i is wall, j is fluid
void ComputeRHEOInterface::store_forces()
{
double *fx = atom->dvector[index_fx];
double *fy = atom->dvector[index_fy];
double *fz = atom->dvector[index_fz];
double minv;
double mass = atom->mass;
double type = atom->type;
double **f = atom->f;
double **fp = atom->fp;
int *mask = atom->mask;
int flag;
// When this is called, fp_store stores the pressure force
// After this method, fp_store instead stores non-pressure forces
// and is also normalized by the particles mass
// If forces are overwritten by a fix, there are no pressure forces
// so just normalize
int ifix = modify->find_fix_by_style("setforce");
if (ifix != -1) {
for (int i = 0; i < atom->nlocal; i++) {
minv = 1.0 / mass[type[i]];
if (mask[i] & modify->fix[ifix]->groupbit) {
fx[i] = f[i][0];
fy[i] = f[i][1];
fz[i] = f[i][2];
fp_store[i][0] = f[i][0] * minv;
fp_store[i][1] = f[i][1] * minv;
fp_store[i][2] = f[i][2] * minv;
} else {
fx[i] = f[i][0] - fp[i][0];
fy[i] = f[i][1] - fp[i][1];
fz[i] = f[i][2] - fp[i][2];
fp_store[i][0] = (f[i][0] - fp_store[i][0]) * 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 {
for (int i = 0; i < atom->nlocal; i++) {
fx[i] = f[i][0] - fp[i][0];
fy[i] = f[i][1] - fp[i][1];
fz[i] = f[i][2] - fp[i][2];
minv = 1.0 / mass[type[i]];
fp_store[i][0] = (f[i][0] - fp_store[i][0]) * 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;
}
}
//Forward comm forces -note only needed here b/c property atom will forward otherwise
// Forward comm forces
comm_forward = 3;
comm_stage = 0;
comm->forward_comm_compute(this);
comm->forward_comm(this);
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */