implement mass keyword to turn off/on change in rmass when diameter is changed

doc/src/fix_adapt.rst

@@ -35,7 +35,7 @@ Syntax
          v_name = variable with name that calculates value of aparam
 
 * zero or more keyword/value pairs may be appended
-* keyword = *scale* or *reset*
+* keyword = *scale* or *reset* or *mass*
 
   .. parsed-literal::
 
@@ -45,6 +45,9 @@ Syntax
        *reset* value = *no* or *yes*
          *no* = values will remain altered at the end of a run
          *yes* = reset altered values to their original values at the end of a run
+       *mass* value = *no* or *yes*
+         *no* = mass is not altered by changes in diameter
+         *yes* = mass is altenred by changes in diameter
 
 Examples
 """"""""
@@ -334,13 +337,16 @@ corresponding attribute for all atoms in the fix group.
 
 If the atom parameter is *diameter* and per-atom density and per-atom
 mass are defined for particles (e.g. :doc:`atom_style granular
-<atom_style>`), then the mass of each particle is also changed when
-the diameter changes. The mass is set from the particle volume for 3d
-systems (density is assumed to stay constant). For 2d, the default is
-for LAMMPS to model particles with a radius attribute as spheres.
-However, if the atom parameter is *diameter/disc*, then the mass is
-set from the particle area (the density is assumed to be in
-mass/distance^2 units).
+<atom_style>`), then the mass of each particle is, by default, also
+changed when the diameter changes. The mass is set from the particle
+volume for 3d systems (density is assumed to stay constant). For 2d,
+the default is for LAMMPS to model particles with a radius attribute
+as spheres. However, if the atom parameter is *diameter/disc*, then the
+mass is set from the particle area (the density is assumed to be in
+mass/distance^2 units). The mass of the particle may also be kept constant
+if the *mass* keyword is set to *no*. This can be useful to account for
+diameter changes that do not involve mass changes, e.g., thermal expansion.
+
 
 For example, these commands would shrink the diameter of all granular
 particles in the "center" group from 1.0 to 0.1 in a linear fashion
@@ -426,4 +432,4 @@ Related commands
 Default
 """""""
 
-The option defaults are scale = no, reset = no.
+The option defaults are scale = no, reset = no, mass = yes.

src/fix_adapt.cpp

@@ -167,6 +167,7 @@ nadapt(0), id_fix_diam(NULL), id_fix_chg(NULL), adapt(NULL)
 
   resetflag = 0;
   scaleflag = 0;
+  massflag = 1;
 
   while (iarg < narg) {
     if (strcmp(arg[iarg],"reset") == 0) {
@@ -181,6 +182,12 @@ nadapt(0), id_fix_diam(NULL), id_fix_chg(NULL), adapt(NULL)
       else if (strcmp(arg[iarg+1],"yes") == 0) scaleflag = 1;
       else error->all(FLERR,"Illegal fix adapt command");
       iarg += 2;
+    } else if (strcmp(arg[iarg],"mass") == 0) {
+      if (iarg+2 > narg)error->all(FLERR,"Illegal fix adapt command");
+      if (strcmp(arg[iarg+1],"no") == 0) massflag = 0;
+      else if (strcmp(arg[iarg+1],"yes") == 0) massflag = 1;
+      else error->all(FLERR,"Illegal fix adapt command");
+      iarg += 2;
     } else error->all(FLERR,"Illegal fix adapt command");
   }
 
@@ -602,12 +609,14 @@ void FixAdapt::change_settings()
 
         for (i = 0; i < nall; i++) {
           if (mask[i] & groupbit) {
+            if (massflag) {
               if (!scaleflag) scale = 0.5*value / radius[i];
-	    if (scaleflag) radius[i] *= scale;
-	    else radius[i] = 0.5*value;
               if (discflag) rmass[i] *= scale*scale;
               else rmass[i] *= scale*scale*scale;
             }
+            if (scaleflag) radius[i] *= scale;
+            else radius[i] = 0.5*value;
+          }
         }
 	  
         if (scaleflag) previous_diam_scale = value;
@@ -704,11 +713,13 @@ void FixAdapt::restore_settings()
 
         for (int i = 0; i < nlocal; i++)
           if (mask[i] & groupbit) {
+            if (massflag) {
               if (!scaleflag) scale = vec[i] / radius[i];
-            radius[i] = vec[i];
               if (discflag) rmass[i] *= scale*scale;
               else rmass[i] *= scale*scale*scale;
             }
+            radius[i] = vec[i];
+          }
       }
       if (chgflag) {
         double *vec = fix_chg->vstore;

src/fix_adapt.h

@@ -44,7 +44,7 @@ class FixAdapt : public Fix {
   void restart(char *);
 
  private:
-  int nadapt,resetflag,scaleflag;
+  int nadapt,resetflag,scaleflag,massflag;
   int anypair, anybond;
   int nlevels_respa;
   char *id_fix_diam,*id_fix_chg;