debugging surface compute

doc/src/fix_rheo.rst

@@ -20,9 +20,10 @@ Syntax
 
        *shift* values = none, turns on velocity shifting
        *thermal* values = none, turns on thermal evolution
-       *surface/detection* values = *sdstyle* *limit*
+       *surface/detection* values = *sdstyle* *limit* *limit/splash*
          *sdstyle* = *coordination* or *divergence*
          *limit* = threshold for surface particles (unitless)
+         *limit/splash* = threshold for splash particles (unitless)
        *interface/reconstruct* values = none, reconstructs interfaces with solid particles
        *rho/sum* values = none, uses the kernel to compute the density of particles
        *density* values = *rho0* (density)

src/RHEO/compute_rheo_surface.cpp

@@ -82,6 +82,8 @@ void ComputeRHEOSurface::init()
   threshold_style = fix_rheo->surface_style;
   threshold_divr = fix_rheo->divr_surface;
   threshold_z = fix_rheo->zmin_surface;
+  threshold_splash = fix_rheo->zmin_splash;
+  interface_flag = fix_rheo->interface_flag;
 
   cutsq = cut * cut;
 
@@ -202,6 +204,7 @@ void ComputeRHEOSurface::compute_peratom()
         rhoj = rho[j];
 
         // Add corrections for walls
+        if (interface_flag) {
           if (fluidi && (!fluidj)) {
             rhoj = compute_interface->correct_rho(j, i);
           } else if ((!fluidi) && fluidj) {
@@ -210,6 +213,7 @@ void ComputeRHEOSurface::compute_peratom()
             rhoi = rho0;
             rhoj = rho0;
           }
+        }
 
         Voli = mass[itype] / rhoi;
         Volj = mass[jtype] / rhoj;
@@ -262,7 +266,7 @@ void ComputeRHEOSurface::compute_peratom()
         if (divr[i] < threshold_divr) {
           status[i] |= STATUS_SURFACE;
           rsurface[i] = 0.0;
-          if (coordination[i] < threshold_z)
+          if (coordination[i] < threshold_splash)
             status[i] |= STATUS_SPLASH;
         }
       }
@@ -272,10 +276,10 @@ void ComputeRHEOSurface::compute_peratom()
       if (mask[i] & groupbit) {
         status[i] |= STATUS_BULK;
         rsurface[i] = cut;
-        if (coordination[i] < threshold_divr) {
+        if (coordination[i] < threshold_z) {
           status[i] |= STATUS_SURFACE;
           rsurface[i] = 0.0;
-          if (coordination[i] < threshold_z)
+          if (coordination[i] < threshold_splash)
             status[i] |= STATUS_SPLASH;
         }
       }

src/RHEO/compute_rheo_surface.h

@@ -41,7 +41,7 @@ class ComputeRHEOSurface : public Compute {
 
  private:
   double cut, cutsq, rho0, threshold_divr;
-  int surface_style, nmax_store, threshold_z;
+  int surface_style, nmax_store, threshold_z, threshold_splash, interface_flag;
   double **B, **gradC;
   int threshold_style, comm_stage;
 

src/RHEO/fix_rheo.cpp

@@ -97,14 +97,16 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
       if (strcmp(arg[iarg + 1], "coordination")) {
         surface_style = COORDINATION;
         zmin_surface = utils::inumeric(FLERR,arg[iarg + 2],false,lmp);
+        zmin_splash = utils::inumeric(FLERR,arg[iarg + 3],false,lmp);
       } else if (strcmp(arg[iarg + 1], "divergence")) {
         surface_style = DIVR;
         divr_surface = utils::numeric(FLERR,arg[iarg + 2],false,lmp);
+        zmin_splash = utils::inumeric(FLERR,arg[iarg + 3],false,lmp);
       } else {
         error->all(FLERR,"Illegal surface/detection option in fix rheo, {}", arg[iarg + 1]);
       }
 
-      iarg += 2;
+      iarg += 3;
     } else if (strcmp(arg[iarg],"interface/reconstruct") == 0) {
       interface_flag = 1;
     } else if (strcmp(arg[iarg],"rho/sum") == 0) {

src/RHEO/fix_rheo.h

@@ -40,7 +40,7 @@ class FixRHEO : public Fix {
 
   // Model parameters
   double h, cut, rho0, csq;
-  int zmin_kernel, zmin_surface;
+  int zmin_kernel, zmin_surface, zmin_splash;
   int kernel_style, surface_style;
   double divr_surface;