From 6de50fbd33a130e8a2ecac07a5bc31ac5467969d Mon Sep 17 00:00:00 2001 From: jtclemm Date: Wed, 14 Jun 2023 13:39:24 -0600 Subject: [PATCH] Adding iso options to fix deform --- doc/src/fix_deform.rst | 46 ++++++++++++--- src/fix_deform.cpp | 129 ++++++++++++++++++++++++++++++++++++----- src/fix_deform.h | 2 + 3 files changed, 155 insertions(+), 22 deletions(-) diff --git a/doc/src/fix_deform.rst b/doc/src/fix_deform.rst index d46f1204f5..7a17f14fed 100644 --- a/doc/src/fix_deform.rst +++ b/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 -evolves according to the equation +specified target value using a linear controller where the box length +: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 ` option. @@ -512,6 +517,25 @@ of the applied strain using the :ref:`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 -only applies to the *pressure* and *pressure/mean* options. +maximum magnitude of the instantaneous strain rate applied in any dimension. +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: diff --git a/src/fix_deform.cpp b/src/fix_deform.cpp index 29c286ed02..e76b34f302 100644 --- a/src/fix_deform.cpp +++ b/src/fix_deform.cpp @@ -65,8 +65,8 @@ irregular(nullptr), set(nullptr) // set defaults - set = new Set[6]; - memset(set,0,6*sizeof(Set)); + set = new Set[7]; + 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; diff --git a/src/fix_deform.h b/src/fix_deform.h index eda97f7c90..e8a4766b12 100644 --- a/src/fix_deform.h +++ b/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(); };