Adding iso options to fix deform

doc/src/fix_deform.rst

@@ -20,7 +20,7 @@ Syntax
 
   .. parsed-literal::
 
-     parameter = *x* or *y* or *z* or *xy* or *xz* or *yz*
+     parameter = *x* or *y* or *z* or *xy* or *xz* or *yz* or *iso*
        *x*, *y*, *z* args = style value(s)
          style = *final* or *delta* or *scale* or *vel* or *erate* or *trate* or *volume* or *wiggle* or *variable*
            *final* values = lo hi
@@ -71,6 +71,12 @@ Syntax
            *variable* values = v_name1 v_name2
              v_name1 = variable with name1 for tilt change as function of time
              v_name2 = variable with name2 for change rate as function of time
+       *iso* = style value
+         style = *volume* or *pressure*
+           *volume* value = none = isotropically adjust system to preserve volume of system
+           *pressure* values = target gain
+             target = target mean pressure (pressure units)
+             gain = proportional gain constant (1/(time * pressure) or 1/time units)
 
 * zero or more keyword/value pairs may be appended
 * keyword = *remap* or *flip* or *units* or *couple* or *vol/balance/p* or *max/rate* or *normalize/pressure*
@@ -88,12 +94,10 @@ Syntax
          box = distances are defined in simulation box units
        *couple* value = *none* or *xyz* or *xy* or *yz* or *xz*
          couple pressure values of various dimensions
-       *vol/balance/p* = *yes* or *no*
+       *vol/balance/p* value = *yes* or *no*
          Modifies the behavior of the *volume* option to try and balance pressures
        *max/rate* value = *rate*
          rate = maximum strain rate for pressure control
-       *normalize/pressure* value = *yes* or *no*
-         determine whether pressure deviation is normalized by target pressure
 
 Examples
 """"""""
@@ -107,6 +111,7 @@ Examples
    fix 1 all deform 1 x pressure 2.0 0.1 normalize/pressure yes max/rate 0.001
    fix 1 all deform 1 x trate 0.1 y volume z volume vol/balance/p yes
    fix 1 all deform 1 x trate 0.1 y pressure/mean 0.0 1.0 z pressure/mean 0.0 1.0
+   fix 1 all deform 1 x trate 0.1 y trate -0.1 overlay/pressure/mean 1.0 0.1
 
 Description
 """""""""""
@@ -318,8 +323,8 @@ mid point.
 
 The *pressure* style adjusts a dimensions's box length to control that
 component of the pressure tensor. This option attempts to maintain a
-specified target value using a linear controller where the box length L
+specified target value using a linear controller where the box length
-evolves according to the equation
+:math:`L` evolves according to the equation
 
 .. parsed-literal::
 
@@ -343,7 +348,7 @@ option.
 
 The *pressure/mean* style is changes a dimension in order to maintain
 a constant mean pressure defined as the trace of the pressure tensor.
-This option is therefore very similar to the *presssure* style with
+This option is therefore very similar to the *pressure* style with
 identical arguments except the current and target pressures refer to the
 mean trace of the pressure tensor. The same options also apply except
 for the :ref:`couple <deform_couple>` option.
@@ -512,6 +517,25 @@ of the applied strain using the :ref:`max/rate <deform_max_rate>` option.
 
 ----------
 
+The *iso* parameter provides an additonal control over the x, y,
+and z box lengths. This parameter can only be used in combination with
+the *x*, *y*, or *z* comamnds: *vel*, *erate*, *trate*, *pressure*, or
+*wiggle*. Note that this parameter will change the overall strain rate in
+the *x*, *y*, or *z* dimensions. This is the meaning of its styles and values.
+
+The *volume* style isotropically scales box lengths to maintain a constant
+box volume in response to deformation from other parameters.
+
+The *pressure* style controls the box volume to maintain the mean pressure
+of the system. This is accomplished by isotropically scaling all box
+lengths :math:`L` by an additional factor of :math:`k (P_t - P_m)` where
+:math:`k` is the proportional gain constant, :math:`P_t` is the target
+pressure, and :math:`P_m` is the current mean pressure (the trace of the
+pressure tensor). This style allows one to control the deviatoric strain
+tensor while maintaining a fixed mean pressure.
+
+----------
+
 All of the tilt styles change the xy, xz, yz tilt factors during a
 simulation.  In LAMMPS, tilt factors (xy,xz,yz) for triclinic boxes
 are normally bounded by half the distance of the parallel box length.
@@ -652,8 +676,12 @@ described below, which will cap the divergence.
 .. _deform_max_rate:
 
 The *max/rate* keyword sets an upper threshold, *rate*, that limits the
-maximum magnitude of the strain rate applied in any dimension. This keyword
+maximum magnitude of the instantaneous strain rate applied in any dimension.
-only applies to the *pressure* and *pressure/mean* options.
+This keyword only applies to the *pressure* and *pressure/mean* options. If
+a pressure-controlled rate is used for both *iso* and either *x*, *y*, or
+*z*, then this threshold will apply separately to each individual controller
+such that the cumulative strain rate on a box dimension may be up to twice
+the value of *rate*.
 
 .. _deform_couple:
 

src/fix_deform.cpp

@@ -65,8 +65,8 @@ irregular(nullptr), set(nullptr)
 
   // set defaults
 
-  set = new Set[6];
+  set = new Set[7];
-  memset(set,0,6*sizeof(Set));
+  memset(set,0,7*sizeof(Set));
 
   // parse arguments
 
@@ -237,6 +237,25 @@ irregular(nullptr), set(nullptr)
           error->all(FLERR,"Illegal fix deform pressure gain, must be positive");
         iarg += 4;
       } else error->all(FLERR,"Illegal fix deform command: {}", arg[iarg+1]);
+    } else if (strcmp(arg[iarg],"iso") == 0) {
+      index = 6;
+      if (strcmp(arg[iarg+1],"volume") == 0) {
+        set[index].style = VOLUME;
+        iarg += 2;
+      } else if (strcmp(arg[iarg+1],"pressure") == 0) {
+        if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "fix deform pressure", error);
+        set[index].style = PRESSURE;
+        if (strstr(arg[iarg+2],"v_") != arg[iarg+2]) {
+          set[index].ptarget = utils::numeric(FLERR,arg[iarg+2],false,lmp);
+        } else {
+          set[index].pstr = utils::strdup(&arg[iarg+2][2]);
+          set[index].pvar_flag = 1;
+        }
+        set[index].pgain = utils::numeric(FLERR,arg[iarg+3],false,lmp);
+        if (set[index].pgain <= 0.0)
+          error->all(FLERR,"Illegal fix deform pressure gain, must be positive");
+        iarg += 4;
+      }
     } else break;
   }
 
@@ -313,6 +332,12 @@ irregular(nullptr), set(nullptr)
     if (set[i].style == NONE) dimflag[i] = 0;
     else dimflag[i] = 1;
 
+  if (set[6].style != NONE) {
+    dimflag[0] = 1;
+    dimflag[1] = 1;
+    dimflag[2] = 1;
+  }
+
   if (dimflag[0]) box_change |= BOX_CHANGE_X;
   if (dimflag[1]) box_change |= BOX_CHANGE_Y;
   if (dimflag[2]) box_change |= BOX_CHANGE_Z;
@@ -324,7 +349,7 @@ irregular(nullptr), set(nullptr)
   // b/c shrink wrap will change box-length
 
   for (int i = 0; i < 3; i++)
-    if (set[i].style && (domain->boundary[i][0] >= 2 || domain->boundary[i][1] >= 2))
+    if ((set[i].style || set[6].style) && (domain->boundary[i][0] >= 2 || domain->boundary[i][1] >= 2))
       error->all(FLERR,"Cannot use fix deform on a shrink-wrapped boundary");
 
   // no tilt deformation on shrink-wrapped 2nd dim
@@ -431,7 +456,7 @@ irregular(nullptr), set(nullptr)
   // set varflag
 
   varflag = 0;
-  for (int i = 0; i < 6; i++) {
+  for (int i = 0; i < 7; i++) {
     if (set[i].style == VARIABLE) varflag = 1;
     if (set[i].pvar_flag) varflag = 1;
   }
@@ -439,7 +464,7 @@ irregular(nullptr), set(nullptr)
   // set pressure_flag
 
   pressure_flag = 0;
-  for (int i = 0; i < 6; i++) {
+  for (int i = 0; i < 7; i++) {
     if (set[i].style == PRESSURE || set[i].style == PMEAN) pressure_flag = 1;
     if (set[i].coupled_flag) pressure_flag = 1;
   }
@@ -452,6 +477,14 @@ irregular(nullptr), set(nullptr)
       if (set[i].style == PMEAN)
         error->all(FLERR, "Cannot use fix deform to assign constant volume and pressure");
 
+  // check conflicts between x,y,z styles and iso
+
+  if (set[6].style)
+    for (int i = 0; i < 3; i++) {
+      if (set[i].style == FINAL || set[i].style == DELTA || set[i].style == SCALE || set[i].style == PMEAN || set[i].style == VARIABLE)
+        error->all(FLERR, "Cannot use fix deform iso parameter with x, y, or z styles other than vel, erate, trate, pressure, and wiggle");
+    }
+
   // check pressure used for max rate and normalize error flag
 
   if (!pressure_flag && max_h_rate != 0)
@@ -470,6 +503,7 @@ irregular(nullptr), set(nullptr)
   set[3].tilt_initial = domain->yz;
   set[4].tilt_initial = domain->xz;
   set[5].tilt_initial = domain->xy;
+  set[6].vol_initial = domain->xprd * domain->yprd * domain->zprd;
 
   // reneighboring only forced if flips can occur due to shape changes
 
@@ -511,7 +545,7 @@ irregular(nullptr), set(nullptr)
 FixDeform::~FixDeform()
 {
   if (set) {
-    for (int i = 0; i < 6; i++) {
+    for (int i = 0; i < 7; i++) {
       delete[] set[i].hstr;
       delete[] set[i].hratestr;
       delete[] set[i].pstr;
@@ -569,7 +603,7 @@ void FixDeform::init()
 
   // check variables for VARIABLE style
 
-  for (int i = 0; i < 6; i++) {
+  for (int i = 0; i < 7; i++) {
     if (set[i].style != VARIABLE) continue;
     set[i].hvar = input->variable->find(set[i].hstr);
     if (set[i].hvar < 0)
@@ -585,7 +619,7 @@ void FixDeform::init()
 
   // check optional variables for PRESSURE or PMEAN style
 
-  for (int i = 0; i < 6; i++) {
+  for (int i = 0; i < 7; i++) {
     if (!set[i].pvar_flag) continue;
     set[i].pvar = input->variable->find(set[i].pstr);
     if (set[i].pvar < 0)
@@ -701,6 +735,8 @@ void FixDeform::init()
     }
   }
 
+  set[6].vol_start = domain->xprd * domain->yprd * domain->zprd;
+
   // if using tilt TRATE, then initial tilt must be non-zero
 
   for (int i = 3; i < 6; i++)
@@ -740,6 +776,7 @@ void FixDeform::init()
   // set domain->h_rate values for use by domain and other fixes/computes
   // initialize all rates to 0.0
   // cannot set here for TRATE,VOLUME,WIGGLE,VARIABLE,PRESSURE since not constant
+  // if iso style is used, these will also not be constant
 
   h_rate = domain->h_rate;
   h_ratelo = domain->h_ratelo;
@@ -867,6 +904,10 @@ void FixDeform::end_of_step()
 
   if (volume_flag) set_volume();
 
+  // apply any final isotropic scalings
+
+  if (set[6].style) set_iso();
+
   // Save pressure/strain rate if required
 
   if (pressure_flag) {
@@ -983,15 +1024,15 @@ void FixDeform::end_of_step()
   // reset global and local box to new size/shape
   // only if deform fix is controlling the dimension
 
-  if (set[0].style) {
+  if (set[0].style || set[6].style) {
     domain->boxlo[0] = set[0].lo_target;
     domain->boxhi[0] = set[0].hi_target;
   }
-  if (set[1].style) {
+  if (set[1].style || set[6].style) {
     domain->boxlo[1] = set[1].lo_target;
     domain->boxhi[1] = set[1].hi_target;
   }
-  if (set[2].style) {
+  if (set[2].style || set[6].style) {
     domain->boxlo[2] = set[2].lo_target;
     domain->boxhi[2] = set[2].hi_target;
   }
@@ -1285,6 +1326,67 @@ void FixDeform::set_volume()
   }
 }
 
+/* ----------------------------------------------------------------------
+   apply isotropic controls
+------------------------------------------------------------------------- */
+
+void FixDeform::set_iso()
+{
+  int i;
+  double scale, shift;
+  double v_rate;
+
+  if (set[6].style == VOLUME) {
+    double v0 = set[6].vol_start;
+    double v = 1.0;
+    for (i = 0; i < 3; i++)
+      v *= (set[i].hi_target - set[i].lo_target);
+
+    scale = std::pow(v0 / v, THIRD);
+    for (i = 0; i < 3; i++) {
+      shift = 0.5 * (set[i].hi_target - set[i].lo_target) * scale;
+      set[i].lo_target = 0.5 * (set[i].lo_start + set[i].hi_start) - shift;
+      set[i].hi_target = 0.5 * (set[i].lo_start + set[i].hi_start) + shift;
+
+      // Recalculate h_rate
+      h_rate[i] = (set[i].hi_target - set[i].lo_target) / (domain->boxhi[i] - domain->boxlo[i]) - 1.0;
+      h_rate[i] /= update->dt;
+    }
+
+  } else if (set[6].style == PRESSURE) {
+
+    // If variable pressure, calculate current target
+    if (set[6].pvar_flag)
+      set[6].ptarget = input->variable->compute_equal(set[6].pvar);
+
+    v_rate = set[6].pgain * (pressure->scalar- set[6].ptarget);
+
+    if (normalize_pressure_flag) {
+      if (set[6].ptarget == 0) {
+        if (max_h_rate == 0) {
+          error->all(FLERR, "Cannot normalize error for zero pressure without defining a max rate");
+        } else v_rate = max_h_rate * v_rate / fabs(v_rate);
+      } else v_rate /= fabs(set[6].ptarget);
+    }
+
+    if (max_h_rate != 0)
+      if (fabs(set[6].ptarget) > max_h_rate)
+        v_rate = max_h_rate * v_rate / fabs(v_rate);
+
+    set[6].cumulative_strain += update->dt * v_rate;
+    scale = (1.0 + set[6].cumulative_strain);
+    for (i = 0; i < 3; i++) {
+      shift = 0.5 * (set[i].hi_target - set[i].lo_target) * scale;
+      set[i].lo_target = 0.5 * (set[i].lo_start + set[i].hi_start) - shift;
+      set[i].hi_target = 0.5 * (set[i].lo_start + set[i].hi_start) + shift;
+
+      // Recalculate h_rate
+      h_rate[i] = (set[i].hi_target - set[i].lo_target) / (domain->boxhi[i] - domain->boxlo[i]) - 1.0;
+      h_rate[i] /= update->dt;
+    }
+  }
+}
+
 /* ----------------------------------------------------------------------
    write Set data to restart file
 ------------------------------------------------------------------------- */
@@ -1294,7 +1396,7 @@ void FixDeform::write_restart(FILE *fp)
   if (comm->me == 0) {
     int size = 6 * sizeof(Set);
     fwrite(&size,sizeof(int),1,fp);
-    fwrite(set,sizeof(Set),6,fp);
+    fwrite(set,sizeof(Set),7,fp);
   }
 }
 
@@ -1306,7 +1408,7 @@ void FixDeform::restart(char *buf)
 {
   int samestyle = 1;
   Set *set_restart = (Set *) buf;
-  for (int i=0; i<6; ++i) {
+  for (int i=0; i<7; ++i) {
     // restore data from initial state
     set[i].lo_initial = set_restart[i].lo_initial;
     set[i].hi_initial = set_restart[i].hi_initial;
@@ -1315,6 +1417,7 @@ void FixDeform::restart(char *buf)
     set[i].saved = set_restart[i].saved;
     set[i].prior_rate = set_restart[i].prior_rate;
     set[i].prior_pressure = set_restart[i].prior_pressure;
+    set[i].cumulative_strain = set_restart[i].cumulative_strain;
     // check if style settings are consistent (should do the whole set?)
     if (set[i].style != set_restart[i].style)
       samestyle = 0;

src/fix_deform.h

@@ -72,6 +72,7 @@ class FixDeform : public Fix {
     double vol_initial, vol_start;
     double ptarget, pgain;
     double prior_pressure, prior_rate;
+    double cumulative_strain;
     int saved;
     int fixed, dynamic1, dynamic2;
     char *hstr, *hratestr, *pstr;
@@ -85,6 +86,7 @@ class FixDeform : public Fix {
   void set_strain();
   void set_pressure();
   void set_volume();
+  void set_iso();
   void couple();
 };