From bad1188c526672b2055caf6be5e01da1875c89b6 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Wed, 22 Mar 2023 17:27:51 -0600 Subject: [PATCH] Reorganizing peratom arrays --- src/RHEO/compute_rheo_grad.cpp | 160 ++++++++++++++++---------------- src/RHEO/compute_rheo_grad.h | 4 +- src/RHEO/fix_rheo.cpp | 2 +- src/RHEO/fix_rheo_thermal.cpp | 49 ++++++---- src/RHEO/fix_rheo_thermal.h | 2 +- src/RHEO/fix_rheo_viscosity.cpp | 44 ++++++--- src/RHEO/fix_rheo_viscosity.h | 6 +- src/RHEO/pair_rheo.cpp | 13 ++- 8 files changed, 162 insertions(+), 118 deletions(-) diff --git a/src/RHEO/compute_rheo_grad.cpp b/src/RHEO/compute_rheo_grad.cpp index e87d39aa53..f963485157 100644 --- a/src/RHEO/compute_rheo_grad.cpp +++ b/src/RHEO/compute_rheo_grad.cpp @@ -23,11 +23,9 @@ #include "force.h" #include "neighbor.h" #include "neigh_list.h" -#include "modify.h" #include "update.h" #include -#include 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,39 +307,33 @@ 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,33 +353,32 @@ 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,22 +393,21 @@ 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; +} \ No newline at end of file diff --git a/src/RHEO/compute_rheo_grad.h b/src/RHEO/compute_rheo_grad.h index 52c5d7c924..220b5813e3 100644 --- a/src/RHEO/compute_rheo_grad.h +++ b/src/RHEO/compute_rheo_grad.h @@ -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; diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index e2819ad320..d142af6483 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -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; diff --git a/src/RHEO/fix_rheo_thermal.cpp b/src/RHEO/fix_rheo_thermal.cpp index 406439f93e..f000bf65cc 100644 --- a/src/RHEO/fix_rheo_thermal.cpp +++ b/src/RHEO/fix_rheo_thermal.cpp @@ -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++) { diff --git a/src/RHEO/fix_rheo_thermal.h b/src/RHEO/fix_rheo_thermal.h index c0b5255caa..e4ed426fb4 100644 --- a/src/RHEO/fix_rheo_thermal.h +++ b/src/RHEO/fix_rheo_thermal.h @@ -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; diff --git a/src/RHEO/fix_rheo_viscosity.cpp b/src/RHEO/fix_rheo_viscosity.cpp index f57410783c..3dfbfc1058 100644 --- a/src/RHEO/fix_rheo_viscosity.cpp +++ b/src/RHEO/fix_rheo_viscosity.cpp @@ -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; diff --git a/src/RHEO/fix_rheo_viscosity.h b/src/RHEO/fix_rheo_viscosity.h index 31c1441e40..d531507a72 100644 --- a/src/RHEO/fix_rheo_viscosity.h +++ b/src/RHEO/fix_rheo_viscosity.h @@ -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; + 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; }; diff --git a/src/RHEO/pair_rheo.cpp b/src/RHEO/pair_rheo.cpp index c796db208a..fd4e525b34 100644 --- a/src/RHEO/pair_rheo.cpp +++ b/src/RHEO/pair_rheo.cpp @@ -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;