From 99e7673e8ea44429a52c3b0c53f68dde001b75a4 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Tue, 21 Feb 2023 19:41:52 -0700 Subject: [PATCH] Starting pair style, various clean ups --- src/RHEO/fix_rheo.cpp | 23 +++++------ src/RHEO/fix_rheo.h | 63 +++++++++++++++-------------- src/RHEO/fix_rheo_thermal.cpp | 70 +++++++++++++++++++++++---------- src/RHEO/fix_rheo_viscosity.cpp | 12 +++--- src/RHEO/pair_rheo.cpp | 58 ++++++++------------------- src/RHEO/pair_rheo.h | 15 ++++--- 6 files changed, 126 insertions(+), 115 deletions(-) diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index 5ff72045f7..e2819ad320 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -66,7 +66,7 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : if (narg < 5) error->all(FLERR,"Insufficient arguments for fix rheo command"); - cut = utils::numeric(FLERR,arg[3],false,lmp); + h = utils::numeric(FLERR,arg[3],false,lmp); if (strcmp(arg[4],"Quintic") == 0) { kernel_style = QUINTIC; } else if (strcmp(arg[4],"CRK0") == 0) { @@ -125,11 +125,14 @@ FixRHEO::~FixRHEO() if (compute_vshift) modify->delete_compute("rheo_vshift"); } -/* ---------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Create necessary internal computes +------------------------------------------------------------------------- */ void FixRHEO::post_constructor() { - compute_kernel = dynamic_cast(modify->add_compute(fmt::format("rheo_kernel all rheo/kernel {} {} {}", kernel_style, zmin_kernel, cut))); + compute_kernel = dynamic_cast(modify->add_compute(fmt::format("rheo_kernel all rheo/kernel {} {} {}", kernel_style, zmin_kernel, h))); fix_store_visc = dynamic_cast(modify->add_fix("rheo_store_visc all STORE/PERATOM 0 1")) fix_store_visc->disable = 1; @@ -141,14 +144,14 @@ void FixRHEO::post_constructor() std::string cmd = "rheo_grad all rheo/grad {} velocity rho viscosity"; if (thermal_flag) cmd += "temperature"; - compute_grad = dynamic_cast(modify->add_compute(fmt::format(cmd, cut))); + compute_grad = dynamic_cast(modify->add_compute(fmt::format(cmd, h))); compute_grad->fix_rheo = this; if (rhosum_flag) - compute_rhosum = dynamic_cast(modify->add_compute(fmt::format("rheo_rhosum all rheo/rho/sum {} {}", cut, zmin_rhosum))); + compute_rhosum = dynamic_cast(modify->add_compute(fmt::format("rheo_rhosum all rheo/rho/sum {} {}", h, zmin_rhosum))); if (shift_flag) - compute_vshift = dynamic_cast(modify->add_compute(fmt::format("rheo_vshift all rheo/vshift {}", cut))); + compute_vshift = dynamic_cast(modify->add_compute(fmt::format("rheo_vshift all rheo/vshift {}", h))); if (surface_flag) { fix_store_surf = dynamic_cast(modify->add_fix("rheo_store_surf all STORE/PERATOM 0 1")) @@ -163,7 +166,7 @@ void FixRHEO::post_constructor() } if (interface_flag) { - compute_interface = dynamic_cast(modify->add_compute(fmt::format("rheo_interface all rheo/interface {}", cut))); + compute_interface = dynamic_cast(modify->add_compute(fmt::format("rheo_interface all rheo/interface {}", h))); fix_store_fp = dynamic_cast(modify->add_fix("rheo_store_fp all STORE/PERATOM 0 3")) f_pressure = fix_store_fp->astore; @@ -197,8 +200,7 @@ void FixRHEO::init() void FixRHEO::setup_pre_force(int /*vflag*/) { - // Check to confirm no accessory fixes are yet defined - // FixRHEO must be the first fix + // Check to confirm accessory fixes do not preceed FixRHEO // Note: these fixes set this flag in setup_pre_force() if (viscosity_fix_defined || pressure_fix_defined || thermal_fix_defined || surface_fix_defined) error->all(FLERR, "Fix RHEO must be defined before all other RHEO fixes"); @@ -210,8 +212,7 @@ void FixRHEO::setup_pre_force(int /*vflag*/) void FixRHEO::setup() { - // Check to confirm all accessory fixes are defined - // Does not ensure fixes correctly cover all atoms (could be a subset group) + // Confirm all accessory fixes are defined, may not cover group all // Note: these fixes set this flag in setup_pre_force() if (!viscosity_fix_defined) error->all(FLERR, "Missing fix rheo/viscosity"); diff --git a/src/RHEO/fix_rheo.h b/src/RHEO/fix_rheo.h index 432d319728..3487dd7273 100644 --- a/src/RHEO/fix_rheo.h +++ b/src/RHEO/fix_rheo.h @@ -32,12 +32,45 @@ class FixRHEO : public Fix { void post_constructor() override; void init() override; void setup_pre_force(int) override; + void setup() override; void pre_force(int) override; void initial_integrate(int) override; void final_integrate() override; void reset_dt() override; + // Model parameters + double h, rho0, csq; + int zmin_kernel, rhosum_zmin; int kernel_style; + enum {QUINTIC, CRK0, CRK1, CRK2}; + + // Non-persistent per-atom arrays + int *surface; + double *conductivity, *viscosity, *pressure; + double **f_pressure; + + // 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 = FFFFFFF0; + int surfacemask = FFFFFF0F; + + // Accessory fixes/computes int thermal_flag; int rhosum_flag; int shift_flag; @@ -47,13 +80,9 @@ class FixRHEO : public Fix { int viscosity_fix_defined; int pressure_fix_defined; int thermal_fix_defined; + int interface_fix_defined; int surface_fix_defined; - // Non-persistent per-atom arrays - int *surface; - double *conductivity, *viscosity, *pressure; - double **f_pressure; - class FixStorePeratom *fix_store_visc; class FixStorePeratom *fix_store_pres; class FixStorePeratom *fix_store_cond; @@ -66,31 +95,7 @@ class FixRHEO : public Fix { class ComputeRHEORhoSum *compute_rhosum; class ComputeRHEOVShift *compute_vshift; - enum {QUINTIC, CRK0, CRK1, CRK2}; - enum {LINEAR, CUBIC, TAITWATER}; - - enum { - // Phase status - STATUS_FLUID = 1 << 0, //Need to turn off other options - STATUS_REACTIVE = 1 << 1, - STATUS_SOLID = 1 << 2, - STATUS_FREEZING = 1 << 3 - - // Temporary status options - reset in preforce - STATUS_SHIFT = 1 << 4, - STATUS_NO_FORCE = 1 << 5, - - // Surface status - STATUS_BULK = 1 << 6, //Need to turn off other options - STATUS_LAYER = 1 << 7, - STATUS_SURFACE = 1 << 8, - STATUS_SPLASH = 1 << 9, - }; - protected: - double cut, rho0, csq; - int zmin_kernel, rhosum_zmin; - double dtv, dtf; }; diff --git a/src/RHEO/fix_rheo_thermal.cpp b/src/RHEO/fix_rheo_thermal.cpp index 37290e3951..406439f93e 100644 --- a/src/RHEO/fix_rheo_thermal.cpp +++ b/src/RHEO/fix_rheo_thermal.cpp @@ -240,9 +240,11 @@ void FixRHEOThermal::post_integrate() Tci = Tc_type[type[i]]); } - if (Ti > Tci) { //Need to untoggle other phase options + if (Ti > Tci) { + status[i] &= phasemask; status[i] |= FixRHEO::STATUS_FLUID; } else if (!(status[i] & FixRHEO::STATUS_SOLID)) + status[i] &= phasemask; status[i] |= FixRHEO::STATUS_FREEZING; } } @@ -250,31 +252,59 @@ void FixRHEOThermal::post_integrate() } } +/* ---------------------------------------------------------------------- + Only need to update non-evolving conductivity styles after atoms exchange +------------------------------------------------------------------------- */ -add post neighbor then update preforce below -/* ---------------------------------------------------------------------- */ +void FixRHEOThermal::post_neighbor() +{ + int i; + int *type = atom->type; + double *conductivity = fix_rheo->conductivity; + 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 (conductivity_style == CONSTANT) { + for (i = 0; i < nall; i++) + if (mask[i] & groupbit) conductivity[i] = kappa; + } else if (conductivity_style == TYPE) { + for (i = 0; i < nall; i++) + if (mask[i] & groupbit) conductivity[i] = kappa_type[type[i]]; + } + } +} + +/* ---------------------------------------------------------------------- + Update (and forward) evolving conductivity styles every timestep +------------------------------------------------------------------------- */ void FixRHEOThermal::pre_force(int /*vflag*/) { - double *conductivity = atom->conductivity; - int *mask = atom->mask; - int nlocal = atom->nlocal; + // So far, none exist + //int i; + //double *conductivity = fix_rheo->conductivity; + //int *mask = atom->mask; + //int nlocal = atom->nlocal; - // Calculate non-persistent quantities before pairstyles - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - conductivity[i] = calc_kappa(i); - } - } + //if (first_flag & nmax < atom->nmax) { + // nmax = atom->nmax; + // fix_rheo->fix_store_cond->grow_arrays(nmax); + //} - if (conductivity_style == CONSTANT) { - return kappa; - } else if (conductivity_style == TYPE) { - int itype = atom->type[i]; - return(kappa_type[itype]); - } else { - error->all(FLERR, "Invalid style"); - } + //if (conductivity_style == TBD) { + // for (i = 0; i < nlocal; i++) { + // if (mask[i] & groupbit) { + // } + // } + //} + + //if (last_flag && comm_forward) comm->forward_comm(this); } /* ---------------------------------------------------------------------- */ diff --git a/src/RHEO/fix_rheo_viscosity.cpp b/src/RHEO/fix_rheo_viscosity.cpp index 8ac32cd1dc..f57410783c 100644 --- a/src/RHEO/fix_rheo_viscosity.cpp +++ b/src/RHEO/fix_rheo_viscosity.cpp @@ -131,8 +131,9 @@ void FixRHEOViscosity::setup_pre_force(int /*vflag*/) pre_force(0); } - -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Only need to update non-evolving viscosity styles after atoms exchange +------------------------------------------------------------------------- */ void FixRHEOViscosity::post_neighbor() { @@ -150,7 +151,6 @@ void FixRHEOViscosity::post_neighbor() fix_rheo->fix_store_visc->grow_arrays(nmax); } - // Update non-evolving viscosity styles only after atoms can exchange if (viscosity_style == CONSTANT) { for (i = 0; i < nall; i++) if (mask[i] & groupbit) viscosity[i] = eta; @@ -161,14 +161,15 @@ void FixRHEOViscosity::post_neighbor() } } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Update (and forward) evolving viscosity styles every timestep +------------------------------------------------------------------------- */ void FixRHEOViscosity::pre_force(int /*vflag*/) { int i, a, b; double tmp, gdot; - int *type = atom->type; double *viscosity = fix_rheo->viscosity; int *mask = atom->mask; double **gradv = compute_grad->gradv; @@ -181,7 +182,6 @@ void FixRHEOViscosity::pre_force(int /*vflag*/) fix_rheo->fix_store_visc->grow_arrays(nmax); } - // Update viscosity styles that evolve in time every timestep if (viscosity_style == POWER) { for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { diff --git a/src/RHEO/pair_rheo.cpp b/src/RHEO/pair_rheo.cpp index b0443933e1..c796db208a 100644 --- a/src/RHEO/pair_rheo.cpp +++ b/src/RHEO/pair_rheo.cpp @@ -394,23 +394,23 @@ void PairRHEO::allocate() void PairRHEO::settings(int narg, char **arg) { - if(narg < 1) error->all(FLERR,"Illegal pair_style command"); + if (narg < 1) error->all(FLERR,"Illegal pair_style command"); int iarg = 0; while (iarg < narg) { if (strcmp(arg[iarg], "rho/damp") == 0) { - if (iarg+1 >= narg) error->all(FLERR,"Illegal pair_style command"); + if (iarg + 1 >= narg) error->all(FLERR,"Illegal pair_style command"); rho_damp_flag = 1; - rho_damp = utils::numeric(FLERR,arg[iarg+1],false,lmp); - iarg ++; + rho_damp = utils::numeric(FLERR,arg[iarg + 1],false,lmp); + iarg++; } else if (strcmp(arg[iarg], "artificial/visc") == 0) { - if (iarg+1 >= narg) error->all(FLERR,"Illegal pair_style command"); + if (iarg + 1 >= narg) error->all(FLERR,"Illegal pair_style command"); artificial_visc_flag = 1; av = utils::numeric(FLERR,arg[iarg+1],false,lmp); - iarg ++; - } else error->all(FLERR,"Illegal pair_style command"); + iarg++; + } else error->all(FLERR,"Illegal pair_style command, {}", arg[iarg]); iarg++; } } @@ -421,8 +421,8 @@ void PairRHEO::settings(int narg, char **arg) void PairRHEO::coeff(int narg, char **arg) { - if (narg != 4) - error->all(FLERR,"Incorrect number of args for pair_style llns coefficients"); + if (narg != 2) + error->all(FLERR,"Incorrect number of args for pair_style rheo coefficients"); if (!allocated) allocate(); @@ -430,17 +430,8 @@ void PairRHEO::coeff(int narg, char **arg) utils::bounds(FLERR,arg[0],1, atom->ntypes, ilo, ihi,error); utils::bounds(FLERR,arg[1],1, atom->ntypes, jlo, jhi,error); - double rho0_one = utils::numeric(FLERR,arg[2],false,lmp); - double c_one = utils::numeric(FLERR,arg[3],false,lmp); - - if (c_one != 1.0) error->warning(FLERR, "Need c = 1 for assumption in compute rheo/solids"); - if (rho0_one != 1.0) error->warning(FLERR, "Need rho0 = 1 for assumption in compute rheo/solids"); - int count = 0; for (int i = ilo; i <= ihi; i++) { - rho0[i] = rho0_one; - csq[i] = c_one*c_one; - for (int j = 0; j <= atom->ntypes; j++) { setflag[i][j] = 1; count++; @@ -448,7 +439,7 @@ void PairRHEO::coeff(int narg, char **arg) } if (count == 0) - error->all(FLERR,"Incorrect args for pair llns coefficients"); + error->all(FLERR,"Incorrect args for pair rheo coefficients"); } /* ---------------------------------------------------------------------- @@ -457,26 +448,21 @@ void PairRHEO::coeff(int narg, char **arg) void PairRHEO::setup() { - int flag; - int ifix = modify->find_fix_by_style("rheo"); - if (ifix == -1) error->all(FLERR, "Using pair RHEO without fix RHEO"); - fix_rheo = ((FixRHEO *) modify->fix[ifix]); + auto fixes = modify->get_fix_by_style("rheo"); + if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use fix rheo/viscosity"); + fix_rheo = dynamic_cast(fixes[0]); compute_kernel = fix_rheo->compute_kernel; compute_grad = fix_rheo->compute_grad; - compute_sinterpolation = fix_rheo->compute_sinterpolation; + compute_interface = fix_rheo->compute_interface; thermal_flag = fix_rheo->thermal_flag; + h = fix_rheo->h; + csq = fix_rheo->csq; + rho0 = fix_rheo->rho0; - - - - h = utils::numeric(FLERR,arg[0],false,lmp); - if (h <= 0.0) error->all(FLERR,"Illegal pair_style command"); hinv = 1.0 / h; - hinv3 = 3.0 * hinv; laplacian_order = -1; - if (comm->ghost_velocity == 0) error->all(FLERR,"Pair RHEO requires ghost atoms store velocity"); @@ -505,13 +491,3 @@ double PairRHEO::init_one(int i, int j) return cut[i][j]; } - -/* ---------------------------------------------------------------------- */ - -double PairRHEO::single(int /*i*/, int /*j*/, int /*itype*/, int /*jtype*/, - double /*rsq*/, double /*factor_coul*/, double /*factor_lj*/, double &fforce) -{ - fforce = 0.0; - - return 0.0; -} diff --git a/src/RHEO/pair_rheo.h b/src/RHEO/pair_rheo.h index 9a9d751d45..026bc16527 100644 --- a/src/RHEO/pair_rheo.h +++ b/src/RHEO/pair_rheo.h @@ -12,9 +12,9 @@ ------------------------------------------------------------------------- */ #ifdef PAIR_CLASS - +// clang-format off PairStyle(rheo,PairRHEO) - +// clang-format on #else #ifndef LMP_PAIR_RHEO_H @@ -33,14 +33,13 @@ class PairRHEO : public Pair { void coeff(int, char **) override; void setup() override; double init_one(int, int) override; - double single(int, int, int, int, double, double, double, double &) override; protected: - int laplacian_order; // From fix RHEO - double h, csq*, rho0*; // From fix RHEO + double h, csq, rho0; // From fix RHEO - double hsq, hinv, rho0, av, rho_damp; + double hsq, hinv, av, rho_damp; + int laplacian_order; int artificial_visc_flag; int rho_damp_flag; int thermal_flag; @@ -49,11 +48,11 @@ class PairRHEO : public Pair { class ComputeRHEOKernel *compute_kernel; class ComputeRHEOGrad *compute_grad; - class ComputeRHEOSolidInterpolation *compute_sinterpolation; + class ComputeRHEOInterface *compute_interface; class FixRHEO *fix_rheo; }; -} +} // namespace LAMMPS_NS #endif #endif