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

Reorganizing peratom arrays

Browse Source
This commit is contained in:
jtclemm
2023-03-22 17:27:51 -06:00
parent 99e7673e8e
commit bad1188c52
8 changed files with 162 additions and 118 deletions
Download Patch File Download Diff File Expand all files Collapse all files

160
src/RHEO/compute_rheo_grad.cpp
View File

@ -23,11 +23,9 @@
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "modify.h"
#include "update.h"
#include <cmath>
#include <string>
using namespace LAMMPS_NS;
enum{COMMGRAD, COMMFIELD};
@ -35,11 +33,11 @@ enum{COMMGRAD, COMMFIELD};
/* ---------------------------------------------------------------------- */
ComputeRHEOGrad::ComputeRHEOGrad(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), compute_interface(nullptr), compute_kernel(nullptr)
Compute(lmp, narg, arg), compute_interface(nullptr), compute_kernel(nullptr),
gradv(nullptr), gradr(nullptr), gradt(nullptr), gradn(nullptr)
{
if (narg < 4) error->all(FLERR,"Illegal compute rheo/grad command");
velocity_flag = temperature_flag = rho_flag = eta_flag = 0;
for (int iarg = 3; iarg < narg; iarg ++) {
if (strcmp(arg[iarg],"velocity") == 0) velocity_flag = 1;
@ -53,70 +51,62 @@ ComputeRHEOGrad::ComputeRHEOGrad(LAMMPS *lmp, int narg, char **arg) :
ncomm_field = 0;
comm_reverse = 0;
std::string fix_cmd = "rheo_grad_property_atom all property/atom"
int dim = domain->dimension;
if (velocity_flag) {
ncomm_grad += dim * dim;
ncomm_field += dim;
comm_reverse += dim * dim;
fix_cmd += " d2_gradv 9"
indexv = atom->add_custom("rheo_grad_v", 1, dim * dim);
gradv = atom->darray[indexv];
}
if (rho_flag) {
ncomm_grad += dim;
ncomm_field += 1;
comm_reverse += dim;
fix_cmd += " d2_gradr 3"
indexr = atom->add_custom("rheo_grad_rho", 1, dim);
gradr = atom->darray[indexr];
}
if (temperature_flag) {
ncomm_grad += dim;
ncomm_field += 1;
comm_reverse += dim;
fix_cmd += " d2_gradt 3"
indext= atom->add_custom("rheo_grad_temp", 1, dim);
gradt = atom->darray[indext];
}
if (eta_flag) {
ncomm_grad += dim;
comm_reverse += dim;
fix_cmd += " d2_gradn 3"
indexn = atom->add_custom("rheo_grad_eta", 1, dim);
gradn = atom->darray[indexn];
}
// Create a custom atom property so it works with compute property/atom
// Do not create grow callback as there's no reason to copy/exchange data
// Manually grow if nmax_old exceeded
comm_forward = ncomm_grad;
modify->add_fix(fix_cmd);
int tmp1, tmp2, index;
if (velocity_flag) {
index = atom->find_custom("gradv", tmp1, tmp2);
gradv = atom->darray[index];
}
if (rho_flag) {
index = atom->find_custom("gradr", tmp1, tmp2);
gradr = atom->darray[index];
}
if (temperature_flag) {
index = atom->find_custom("gradt", tmp1, tmp2);
gradt = atom->darray[index];
}
if (eta_flag) {
index = atom->find_custom("gradn", tmp1, tmp2);
gradn = atom->darray[index];
}
nmax_old = 0;
grow_arrays(atom->nmax);
}
/* ---------------------------------------------------------------------- */
ComputeRHEOGrad::~ComputeRHEOGrad()
{
modify->delete_fix("rheo_grad_property_atom");
int dim = domain->dimension;
if (velocity_flag)
atom->remove_custom(indexv, 1, dim * dim);
if (rho_flag)
atom->remove_custom(indexr, 1, dim);
if (temperature_flag)
atom->remove_custom(indext, 1, dim);
if (eta_flag)
atom->remove_custom(indexn, 1, dim);
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOGrad::init()
@ -169,6 +159,8 @@ void ComputeRHEOGrad::compute_peratom()
firstneigh = list->firstneigh;
// initialize arrays
if (nmax > nmax_old) grow_arrays(nmax);
for (i = 0; i < nmax; i++) {
if (velocity_flag) {
for (k = 0; k < dim * dim; k++)
@ -315,41 +307,35 @@ int ComputeRHEOGrad::pack_forward_comm(int n, int *list, double *buf,
j = list[i];
if (comm_stage == COMMGRAD) {
if (velocity_flag){
if (velocity_flag)
for (k = 0; k < dim * dim; k++)
buf[m++] = gradv[j][k];
}
if (rho_flag) {
if (rho_flag)
for (k = 0; k < dim; k++)
buf[m++] = gradr[j][k];
}
if (temperature_flag) {
if (temperature_flag)
for (k = 0; k < dim; k++)
buf[m++] = gradt[j][k];
}
if (eta_flag){
if (eta_flag)
for (k = 0; k < dim; k++)
buf[m++] = gradn[j][k];
}
} else if (comm_stage == COMMFIELD) {
if (velocity_flag) {
if (velocity_flag)
for (k = 0; k < dim; k++)
buf[m++] = v[j][k];
}
if (rho_flag) {
if (rho_flag)
buf[m++] = rho[j];
}
if (temperature_flag) {
if (temperature_flag)
buf[m++] = temperature[j];
}
}
}
return m;
}
@ -367,35 +353,34 @@ void ComputeRHEOGrad::unpack_forward_comm(int n, int first, double *buf)
last = first + n;
for (i = first; i < last; i++) {
if (comm_stage == COMMGRAD) {
if (velocity_flag) {
if (velocity_flag)
for (k = 0; k < dim * dim; k++)
gradv[i][k] = buf[m++];
}
if (rho_flag) {
if (rho_flag)
for (k = 0; k < dim; k++)
gradr[i][k] = buf[m++];
}
if (temperature_flag) {
if (temperature_flag)
for (k = 0; k < dim; k++)
gradt[i][k] = buf[m++];
}
if (eta_flag) {
if (eta_flag)
for (k = 0; k < dim; k++)
gradn[i][k] = buf[m++];
}
} else if (comm_stage == COMMFIELD) {
if (velocity_flag) {
if (velocity_flag)
for (k = 0; k < dim; k++)
v[i][k] = buf[m++];
}
if (rho_flag) {
if (rho_flag)
rho[i] = buf[m++];
}
if (temperature_flag) {
if (temperature_flag)
temperature[i] = buf[m++];
}
}
}
}
/* ---------------------------------------------------------------------- */
@ -408,23 +393,22 @@ int ComputeRHEOGrad::pack_reverse_comm(int n, int first, double *buf)
m = 0;
last = first + n;
for (i = first; i < last; i++) {
if (velocity_flag) {
if (velocity_flag)
for (k = 0; k < dim * dim; k++)
buf[m++] = gradv[i][k];
}
if (rho_flag) {
if (rho_flag)
for (k = 0; k < dim; k++)
buf[m++] = gradr[i][k];
}
if (temperature_flag) {
if (temperature_flag)
for (k = 0; k < dim; k++)
buf[m++] = gradt[i][k];
}
if (eta_flag) {
if (eta_flag)
for (k = 0; k < dim; k++)
buf[m++] = gradn[i][k];
}
}
return m;
}
@ -438,21 +422,39 @@ void ComputeRHEOGrad::unpack_reverse_comm(int n, int *list, double *buf)
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
if (velocity_flag) {
if (velocity_flag)
for (k = 0; k < dim * dim; k++)
gradv[j][k] += buf[m++];
}
if (rho_flag) {
if (rho_flag)
for (k = 0; k < dim; k++)
gradr[j][k] += buf[m++];
}
if (temperature_flag) {
if (temperature_flag)
for (k = 0; k < dim; k++)
gradt[j][k] += buf[m++];
}
if (eta_flag) {
if (eta_flag)
for (k = 0; k < dim; k++)
gradn[j][k] += buf[m++];
}
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOGrad::grow_arrays(int nmax)
{
int dim = domain->dimension;
if (velocity_flag)
memory->grow(gradv, nmax, dim * dim, "atom:rheo_grad_v");
if (rho_flag)
memory->grow(gradr, nmax, dim, "atom:rheo_grad_rho");
if (temperature_flag)
memory->grow(gradt, nmax, dim, "atom:rheo_grad_temp");
if (eta_flag)
memory->grow(gradn, nmax, dim, "atom:rheo_grad_eta");
nmax_old = nmax;
}

4
src/RHEO/compute_rheo_grad.h
View File

@ -44,8 +44,8 @@ class ComputeRHEOGrad : public Compute {
int stage;
private:
int comm_stage;
int ncomm_grad, ncomm_field;
int comm_stage, ncomm_grad, ncomm_field, nmax_old;
int indexv, indexr, indext, indexn;
double cut, cutsq, rho0;
class NeighList *list;

2
src/RHEO/fix_rheo.cpp
View File

@ -313,7 +313,7 @@ void FixRHEO::initial_integrate(int /*vflag*/)
// Shifting atoms
if (shift_flag) {
compute_vshift->correct_surfaces();
compute_vshift->correct_surfaces(); // COuld this be moved to preforce after the surface fix runs?
for (i = 0; i < nlocal; i++) {
if (!(status[i] & STATUS_SHIFT)) continue;

49
src/RHEO/fix_rheo_thermal.cpp
View File

@ -32,7 +32,8 @@ enum {NONE, CONSTANT, TYPE};
/* ---------------------------------------------------------------------- */
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), Tc_type(nullptr), kappa_type(nullptr), cv_type(nullptr),
conductivity(nullptr)
{
if (narg < 4) error->all(FLERR,"Illegal fix command");
@ -40,8 +41,8 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
cv_style = NONE;
conductivity_style = NONE;
comm_forward = 0;
nmax = atom->nmax;
comm_forward = 1;
nmax_old = 0;
int ntypes = atom->ntypes;
int iarg = 3;
@ -123,6 +124,11 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
FixRHEOThermal::~FixRHEOThermal()
{
// Remove custom property if it exists
int tmp1, tmp2, index;
index = atom->find_custom("rheo_conductivity", tmp1, tmp2);
if (index != -1) atom->remove_custom(index_cond, 1, 0);
memory->destroy(cv_type);
memory->destroy(Tc_type);
memory->destroy(kappa_type);
@ -164,8 +170,7 @@ void FixRHEOThermal::setup_pre_force(int /*vflag*/)
fix_rheo->thermal_fix_defined = 1;
// Identify whether this is the first/last instance of fix thermal
// First will handle growing arrays
// Last will handle communication
// First will grow arrays, last will communicate
first_flag = 0
last_flag = 0;
@ -179,6 +184,18 @@ void FixRHEOThermal::setup_pre_force(int /*vflag*/)
if (i == 0) first_flag = 1;
if ((i + 1) == fixlist.size()) last_flag = 1;
// Create conductivity array if it doesn't already exist
// Create a custom atom property so it works with compute property/atom
// Do not create grow callback as there's no reason to copy/exchange data
// Manually grow if nmax_old exceeded
int tmp1, tmp2;
index_cond = atom->find_custom("rheo_conductivity", tmp1, tmp2);
if (index_cond == -1) {
index_cond = atom->add_custom("rheo_conductivity", 1, 0);
nmax_old = atom->nmax;
}
post_neighbor();
pre_force(0);
}
@ -230,7 +247,7 @@ void FixRHEOThermal::post_integrate()
if (status[i] == FixRHEO::FLUID_NO_FORCE) continue;
cvi = calc_cv(i);
temperature[i] += dtf*heat[i]/cvi;
temperature[i] += dtf * heat[i] / cvi;
if (Tc_style != NONE) {
Ti = temperature[i];
@ -260,15 +277,14 @@ void FixRHEOThermal::post_neighbor()
{
int i;
int *type = atom->type;
double *conductivity = fix_rheo->conductivity;
double *conductivity = atom->dvector[index_cond];
int *mask = atom->mask;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
if (first_flag & nmax < atom->nmax) {
nmax = atom->nmax;
fix_rheo->fix_store_cond->grow_arrays(nmax);
}
if (first_flag && (nmax_old < atom->nmax))
memory->grow(conductivity, atom->nmax, "atom:rheo_conductivity");
nmax_old = atom->nmax;
if (conductivity_style == CONSTANT) {
for (i = 0; i < nall; i++)
@ -286,16 +302,15 @@ void FixRHEOThermal::post_neighbor()
void FixRHEOThermal::pre_force(int /*vflag*/)
{
// So far, none exist
// Not needed yet, when needed add (un)pack_forward_comm() methods
//int i;
//double *conductivity = fix_rheo->conductivity;
//double *conductivity = atom->dvector[index_cond];
//int *mask = atom->mask;
//int nlocal = atom->nlocal;
//if (first_flag & nmax < atom->nmax) {
// nmax = atom->nmax;
// fix_rheo->fix_store_cond->grow_arrays(nmax);
//}
//if (first_flag && (nmax_old < atom->nmax))
// memory->grow(conductivity, atom->nmax, "atom:rheo_conductivity");
//nmax_old = atom->nmax;
//if (conductivity_style == TBD) {
// for (i = 0; i < nlocal; i++) {

2
src/RHEO/fix_rheo_thermal.h
View File

@ -49,7 +49,7 @@ class FixRHEOThermal : public Fix {
int cv_style;
int conductivity_style;
int first_flag, last_flag;
int nmax;
int nmax_old, index_cond;
class FixRHEO *fix_rheo;

44
src/RHEO/fix_rheo_viscosity.cpp
View File

@ -39,7 +39,7 @@ FixRHEOViscosity::FixRHEOViscosity(LAMMPS *lmp, int narg, char **arg) :
viscosity_style = NONE;
comm_forward = 0;
nmax = atom->nmax;
nmax_old = 0;
int ntypes = atom->ntypes;
int iarg = 3;
@ -81,6 +81,11 @@ FixRHEOViscosity::FixRHEOViscosity(LAMMPS *lmp, int narg, char **arg) :
FixRHEOViscosity::~FixRHEOViscosity()
{
// Remove custom property if it exists
int tmp1, tmp2, index;
index = atom->find_custom("rheo_viscosity", tmp1, tmp2);
if (index != -1) atom->remove_custom(index_visc, 1, 0);
memory->destroy(eta_type);
}
@ -112,8 +117,7 @@ void FixRHEOViscosity::setup_pre_force(int /*vflag*/)
fix_rheo->viscosity_fix_defined = 1;
// Identify whether this is the first/last instance of fix viscosity
// First will handle growing arrays
// Last will handle communication
// First will grow arrays, last will communicate
first_flag = 0
last_flag = 0;
@ -127,6 +131,18 @@ void FixRHEOViscosity::setup_pre_force(int /*vflag*/)
if (i == 0) first_flag = 1;
if ((i + 1) == fixlist.size()) last_flag = 1;
// Create viscosity array if it doesn't already exist
// Create a custom atom property so it works with compute property/atom
// Do not create grow callback as there's no reason to copy/exchange data
// Manually grow if nmax_old exceeded
int tmp1, tmp2;
index_visc = atom->find_custom("rheo_viscosity", tmp1, tmp2);
if (index_cond == -1) {
index_visc = atom->add_custom("rheo_viscosity", 1, 0);
nmax_old = atom->nmax;
}
post_neighbor();
pre_force(0);
}
@ -140,16 +156,15 @@ void FixRHEOViscosity::post_neighbor()
int i;
int *type = atom->type;
double *viscosity = fix_rheo->viscosity;
double *viscosity = atom->dvector[index_visc];
int *mask = atom->mask;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
if (first_flag & nmax < atom->nmax) {
nmax = atom->nmax;
fix_rheo->fix_store_visc->grow_arrays(nmax);
}
if (first_flag && (nmax_old < atom->nmax))
memory->grow(viscosity, atom->nmax, "atom:rheo_viscosity");
nmax_old = atom->nmax;
if (viscosity_style == CONSTANT) {
for (i = 0; i < nall; i++)
@ -170,17 +185,16 @@ void FixRHEOViscosity::pre_force(int /*vflag*/)
int i, a, b;
double tmp, gdot;
double *viscosity = fix_rheo->viscosity;
double *viscosity = atom->dvector[index_visc];
int *mask = atom->mask;
double **gradv = compute_grad->gradv;
int nlocal = atom->nlocal;
int dim = domain->dimension;
if (first_flag & nmax < atom->nmax) {
nmax = atom->nmax;
fix_rheo->fix_store_visc->grow_arrays(nmax);
}
if (first_flag && (nmax_old < atom->nmax))
memory->grow(viscosity, atom->nmax, "atom:rheo_viscosity");
nmax_old = atom->nmax;
if (viscosity_style == POWER) {
for (i = 0; i < nlocal; i++) {
@ -213,7 +227,7 @@ int FixRHEOViscosity::pack_forward_comm(int n, int *list, double *buf,
int /*pbc_flag*/, int * /*pbc*/)
{
int i,j,k,m;
double *viscosity = fix_rheo->viscosity;
double *viscosity = atom->dvector[index_visc];
m = 0;
for (i = 0; i < n; i++) {
@ -228,7 +242,7 @@ int FixRHEOViscosity::pack_forward_comm(int n, int *list, double *buf,
void FixRHEOViscosity::unpack_forward_comm(int n, int first, double *buf)
{
int i, k, m, last;
double *viscosity = fix_rheo->viscosity;
double *viscosity = atom->dvector[index_visc];
m = 0;
last = first + n;

4
src/RHEO/fix_rheo_viscosity.h
View File

@ -35,12 +35,14 @@ class FixRHEOViscosity : public Fix {
void pre_force(int) override;
int pack_forward_comm(int, int *, double *, int, int *) override;
void unpack_forward_comm(int, int, double *) override;
private:
double *eta_type, eta;
double npow, K, gd0, tau0;
int viscosity_style;
int first_flag, last_flag;
int nmax;
int nmax_old, index_visc;
class FixRHEO *fix_rheo;
class ComputeRHEOGrad *compute_grad;
};

13
src/RHEO/pair_rheo.cpp
View File

@ -92,7 +92,6 @@ void PairRHEO::compute(int eflag, int vflag)
double **gradv = compute_grad->gradv;
double **gradt = compute_grad->gradt;
double **gradr = compute_grad->gradr;
double **gradn = compute_grad->gradn;
double **v = atom->v;
double **x = atom->x;
double **f = atom->f;
@ -109,6 +108,18 @@ void PairRHEO::compute(int eflag, int vflag)
int *type = atom->type;
int *phase = atom->phase;
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];
double *conductivity;
if (thermal_flag) {
int index_cond = atom->find_custom("rheo_conductivity", tmp1, tmp2);
if (index_cond == -1) error->all(FLERR, "Cannot find rheo conductivity");
conductivity = atom->dvector[index_cond];
}
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
int dim = domain->dimension;