Revamping atom data storage in fixes

src/RHEO/fix_rheo.cpp

@@ -21,6 +21,7 @@
 #include "compute_rheo_vshift.h"
 #include "domain.h"
 #include "error.h"
+#include "fix_store_peratom.h"
 #include "force.h"
 #include "modify.h"
 #include "update.h"
@@ -33,15 +34,30 @@ using namespace FixConst;
 
 FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg), compute_grad(nullptr), compute_kernel(nullptr),
-  compute_interface(nullptr), compute_rhosum(nullptr), compute_vshift(nullptr)
+  compute_interface(nullptr), compute_rhosum(nullptr), compute_vshift(nullptr),
+  fix_store_visc(nullptr), fix_store_pres(nullptr), fix_store_cond(nullptr),
+  fix_store_surf(nullptr), fix_store_fp(nullptr), surface(nullptr), conductivity(nullptr),
+  viscosity(nullptr), pressure(nullptr), f_pressure(nullptr)
 {
+  time_integrate = 1;
+
+  viscosity_fix_defined = 0;
+  pressure_fix_defined = 0;
+  thermal_fix_defined = 0;
+  surface_fix_defined = 0;
+
   thermal_flag = 0;
   rhosum_flag = 0;
   shift_flag = 0;
-  solid_flag = 0;
+  interface_flag = 0;
+  surface_flag = 0;
+
   rho0 = 1.0;
   csq = 1.0;
 
+  if (igroup != 0)
+    error->all(FLERR,"fix rheo command requires group all");
+
   if (atom->rho_flag != 1)
     error->all(FLERR,"fix rheo command requires atom_style with density");
   if (atom->status_flag != 1)
@@ -68,6 +84,10 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
       shift_flag = 1;
     } else if (strcmp(arg[iarg],"thermal") == 0) {
       thermal_flag = 1;
+    } else if (strcmp(arg[iarg],"surface/detection") == 0) {
+      surface_flag = 1;
+    } else if (strcmp(arg[iarg],"interface/reconstruction") == 0) {
+      interface_flag = 1;
     } else if (strcmp(arg[iarg],"rhosum") == 0) {
       rhosum_flag = 1;
       if(iarg + 1 >= narg) error->all(FLERR,"Illegal rhosum option in fix rheo");
@@ -86,17 +106,18 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
     }
     iarg += 1;
   }
-
-  time_integrate = 1;
-  thermal_fix_defined = 0;
-  viscosity_fix_defined = 0;
-  pressure_fix_defined = 0;
 }
 
 /* ---------------------------------------------------------------------- */
 
 FixRHEO::~FixRHEO()
 {
+  if (fix_store_visc) modify->delete_fix("rheo_store_visc");
+  if (fix_store_pres) modify->delete_fix("rheo_store_pres");
+  if (fix_store_surf) modify->delete_fix("rheo_store_surf");
+  if (fix_store_cond) modify->delete_fix("rheo_store_cond");
+  if (fix_store_fp) modify->delete_fix("rheo_store_fp");
+
   if (compute_kernel) modify->delete_compute("rheo_kernel");
   if (compute_grad) modify->delete_compute("rheo_grad");
   if (compute_interface) modify->delete_compute("rheo_interface");
@@ -110,13 +131,18 @@ void FixRHEO::post_constructor()
 {
   compute_kernel = dynamic_cast<ComputeRHEOKernel *>(modify->add_compute(fmt::format("rheo_kernel all rheo/kernel {} {} {}", kernel_style, zmin_kernel, cut)));
 
-  if (thermal_flag)
+  fix_store_visc = dynamic_cast<FixStorePeratom *>(modify->add_fix("rheo_store_visc all STORE/PERATOM 0 1"))
-    compute_grad = dynamic_cast<ComputeRHEOGrad *>(modify->add_compute(fmt::format("rheo_grad all rheo/grad {}  velocity rho viscosity temprature", cut)));
+  fix_store_visc->disable = 1;
-  else
+  viscosity = fix_store_visc->vstore;
-    compute_grad = dynamic_cast<ComputeRHEOGrad *>(modify->add_compute(fmt::format("rheo_grad all rheo/grad {}  velocity rho viscosity", cut)));
+  fix_store_pres = dynamic_cast<FixStorePeratom *>(modify->add_fix("rheo_store_pres all STORE/PERATOM 0 1"))
-  compute_grad->fix_rheo = this;
+  pressure = fix_store_pres->vstore;
+  fix_store_pres->disable = 1;
 
-  compute_interface = dynamic_cast<ComputeRHEOInterface *>(modify->add_compute(fmt::format("rheo_interface all rheo/interface {}", cut)));
+
+  std::string cmd = "rheo_grad all rheo/grad {} velocity rho viscosity";
+  if (thermal_flag) cmd += "temperature";
+  compute_grad = dynamic_cast<ComputeRHEOGrad *>(modify->add_compute(fmt::format(cmd, cut)));
+  compute_grad->fix_rheo = this;
 
   if (rhosum_flag)
     compute_rhosum = dynamic_cast<ComputeRHEORhoSum *>(modify->add_compute(fmt::format("rheo_rhosum all rheo/rho/sum {} {}", cut, zmin_rhosum)));
@@ -124,11 +150,25 @@ void FixRHEO::post_constructor()
   if (shift_flag)
     compute_vshift = dynamic_cast<ComputeRHEOVShift *>(modify->add_compute(fmt::format("rheo_vshift all rheo/vshift {}", cut)));
 
+  if (surface_flag) {
+    fix_store_surf = dynamic_cast<FixStorePeratom *>(modify->add_fix("rheo_store_surf all STORE/PERATOM 0 1"))
+    surface = fix_store_surf->vstore;
+    fix_store_surf->disable = 1;
+  }
 
-  //todo here
+  if (thermal_flag) {
-  //allocate memory for viscosity, pressure, etc
+    fix_store_cond = dynamic_cast<FixStorePeratom *>(modify->add_fix("rheo_store_cond all STORE/PERATOM 0 1"))
-  //don't want to save to restart/datafiles (could disable fix store/state)
+    conductivity = fix_store_cond->vstore;
-  //but do want it available for dupm files
+    fix_store_cond->disable = 1;
+  }
+
+  if (interface_flag) {
+    compute_interface = dynamic_cast<ComputeRHEOInterface *>(modify->add_compute(fmt::format("rheo_interface all rheo/interface {}", cut)));
+
+    fix_store_fp = dynamic_cast<FixStorePeratom *>(modify->add_fix("rheo_store_fp all STORE/PERATOM 0 3"))
+    f_pressure = fix_store_fp->astore;
+    fix_store_fp->disable = 1;
+  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -158,19 +198,23 @@ void FixRHEO::init()
 void FixRHEO::setup_pre_force(int /*vflag*/)
 {
   // Check to confirm all accessory fixes are defined
-  if(!thermal_fix_defined && thermal_flag)
-    error->all(FLERR, "Missing fix rheo/thermal");
-
   if (!viscosity_fix_defined)
     error->all(FLERR, "Missing fix rheo/viscosity");
 
   if (!pressure_fix_defined)
     error->all(FLERR, "Missing fix rheo/pressure");
 
+  if(!thermal_fix_defined && thermal_flag)
+    error->all(FLERR, "Missing fix rheo/thermal");
+
+  if(!surface_fix_defined && surface_flag)
+    error->all(FLERR, "Missing fix rheo/surface");
+
   // Reset to zero for next run
   thermal_fix_defined = 0;
   viscosity_fix_defined = 0;
   pressure_fix_defined = 0;
+  surface_fix_defined = 0;
 
   pre_force(0);
 }

src/RHEO/fix_rheo.h

@@ -41,17 +41,25 @@ class FixRHEO : public Fix {
   int thermal_flag;
   int rhosum_flag;
   int shift_flag;
-  int solid_flag;
+  int interface_flag;
+  int surface_flag;
 
-  int thermal_fix_defined;
   int viscosity_fix_defined;
   int pressure_fix_defined;
+  int thermal_fix_defined;
+  int surface_fix_defined;
 
-  // Non-persistent per-atom arrays are initialized here
+  // 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;
+  class FixStorePeratom *fix_store_surf;
+  class FixStorePeratom *fix_store_fp;
+
   class ComputeRHEOGrad *compute_grad;
   class ComputeRHEOKernel *compute_kernel;
   class ComputeRHEOInterface *compute_interface;

src/RHEO/fix_rheo_viscosity.cpp

@@ -39,6 +39,7 @@ FixRHEOViscosity::FixRHEOViscosity(LAMMPS *lmp, int narg, char **arg) :
   viscosity_style = NONE;
 
   comm_forward = 1;
+  nmax = atom->nmax;
 
   int ntypes = atom->ntypes;
   int iarg = 3;
@@ -107,8 +108,10 @@ void FixRHEOViscosity::init()
 
 void FixRHEOViscosity::setup_pre_force(int /*vflag*/)
 {
-  // Identify whether this is the last instance of fix viscosity
+  // Identify whether this is the first/last instance of fix viscosity
-  // Will handle communication
+  // First will handle growing arrays
+  // Last will handle communication
+  first_flag = 0
   last_flag = 0;
 
   int i = 0;
@@ -118,6 +121,7 @@ void FixRHEOViscosity::setup_pre_force(int /*vflag*/)
     i++;
   }
 
+  if (i == 0) first_flag = 1;
   if ((i + 1) == fixlist.size()) last_flag = 1;
 
   pre_force(0);
@@ -138,6 +142,11 @@ void FixRHEOViscosity::pre_force(int /*vflag*/)
   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);
+  }
+
   for (i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) {
       if (viscosity_style == CONSTANT) {

src/RHEO/fix_rheo_viscosity.h

@@ -38,7 +38,8 @@ class FixRHEOViscosity : public Fix {
   double *eta_type, eta;
   double  npow, K, gd0, tau0;
   int viscosity_style;
-  int last_flag;
+  int first_flag, last_flag;
+  int nmax;
   class FixRHEO *fix_rheo;
   class ComputeRHEOGrad *compute_grad;
 };